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Thesis  directed  by  Dr.  Scott  E.  Palo 

Due  to  the  volume  and  power  limitations  of  a  small  satellite,  careful  considera¬ 
tion  must  be  taken  while  designing  an  attitude  control  system  for  3-axis  stabilization. 
Placing  redundancy  in  the  system  proves  difficult  and  utilizing  power  hungry,  high 
accuracy,  active  actuators  is  not  a  viable  option.  Thus,  it  is  customary  to  find  depend¬ 
able,  passive  actuators  used  in  conjunction  with  small  scale  active  control  components. 
This  document  describes  the  application  of  Elastic  Memory  Composite  materials  in  the 
construction  of  a  flexible  spacecraft  appendage,  such  as  a  gravity  gradient  boom.  As¬ 
sumed  modes  methods  are  used  with  Finite  Element  Modeling  information  to  obtain 
the  equations  of  motion  for  the  system  while  assuming  free-free  boundary  conditions. 
A  discussion  is  provided  to  illustrate  how  cantilever  mode  shapes  are  not  always  the 
best  assumption  when  modeling  small  flexible  spacecraft.  A  key  point  of  interest  is  first 
resonant  modes  may  be  needed  in  the  system  design  plant  in  spite  of  these  modes  be¬ 
ing  greater  than  one  order  of  magnitude  in  frequency  when  compared  to  the  crossover 
frequency  of  the  controller.  LQG/LTR  optimal  control  techniques  are  implemented 
to  compute  attitude  control  gains  while  controller  robustness  considerations  determine 
appropriate  reduced  order  controllers  and  which  flexible  modes  to  include  in  the  de¬ 
sign  model.  Key  satellite  designer  concerns  in  the  areas  of  computer  processor  sizing, 
material  uncertainty  impacts  on  the  system  model,  and  system  performance  variations 
resulting  from  appendage  length  modifications  are  addressed. 
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Chapter  1 


Introduction 


1.1  Small  Satellites 

The  world  entered  a  new  stage  in  its  technological  advancement  on  October  4th, 
1957.  On  this  day,  Russians  launched  the  first  man  made  satellite  into  orbit.  Sputnik 
was  shaped  like  a  basketball  with  a  diameter  of  2  meters  and  weighed  only  84  kg 
(see  Figure  1.1).  A  month  later,  Sputnik  II,  weighing  a  total  of  511  kg,  placed  a 
dog  named  Laika  in  an  orbit  about  the  Earth.  The  mission  was  declared  a  failure 
when  the  satellite  experienced  difficulty  separating  from  its  booster  and  overheated[123]. 
Seventeen  months  after  the  launch  of  Sputnik  I,  the  United  States  launched  Explorer 
1.  This  light  satellite,  weighing  only  14  kg,  measured  charged  particles  in  the  upper 
atmosphere  and  discovered  the  Van  Allen  radiation  belts  [5]. 

These  satellites  were  small  for  two  reasons.  First,  the  launch  capability  of  existing 
rockets  had  only  been  evolving  for  a  few  decades  prior  to  the  launch  of  Sputnik.  In  fact, 
American  scientists  were  surprised  with  how  well  the  Russians  could  put  heavy  payloads 
(1,500  kg)  into  orbit[186].  Program  risk  is  the  other  reason  for  the  use  of  small  satellites. 
The  first  launch  attempt  by  the  United  States  failed  on  December  6,  1957  when  a  Navy 
Vanguard  rocket  exploded  before  leaving  the  launch  pad.  The  failure  may  have  been 
catastrophic  to  the  American’s  race  for  space  had  they  dedicated  a  majority  of  their 
program  funds  to  launching  a  larger  version  of  Explorer  1. 

Since  the  end  of  the  1950s,  satellites  have  grown  in  mass  and  expense  as  larger 
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Figure  1.1:  Sputnik  I 


missions  were  attempted.  Businesses  emerged  to  take  advantage  of  the  remote  sensing 
and  communication  capabilities  orbiting  satellites  offered.  These  businesses  could  offset 
the  added  risk  of  launching  larger  satellites  by  adding  advanced  instruments  capable  of 
producing  better  results  than  cheaper,  smaller  scale  versions  of  the  payloads.  However, 
advanced  instruments  consume  additional  power  and  the  battery  and  solar  array  mass 
increased  accordingly.  Primary  structures,  which  carry  most  of  the  loads  of  the  satellite, 
also  grew  in  mass  and  complexity  to  support  the  additional  components.  The  Milstar 
Satellite  Communications  System  (shown  in  Figure  1.2)  is  an  example  of  a  large  satel¬ 
lite  currently  in  use.  This  system  provides  secure,  worldwide  communications  to  high 
priority  military  users  to  meet  essential  wartime  requirements.  Each  Lockheed  Martin 
designed  satellite  weighs  4,536  kg  and  costs  approximately  $800  million[47]. 

A  45  year  trend  in  small  satellites  is  shown  in  Figure  1.3.  This  figure  illustrates 
the  number  of  successful  launches  of  small  satellites  with  the  mini,  micro,  and  nano 
classifications  (see  Table  1.1).  Data  is  compiled  from  various  sources  presented  by 
Surrey  Satellite  Technology  Ltd.  (SSTL)  and  does  not  include  nonfunctioning  satellites 
(e.g.  orbital  debris  radar  calibration  spheres  released  by  the  Shuttle  intended  to  provide 
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Figure  1.2:  Milstar  Satellite  Communications  System 


calibration  for  radar  echoes). 

Figure  1.3  begins  with  two  Sputnik  satellites  launched  in  1957.  The  following 
decade  comprises  a  combination  of  US  and  Russian  satellite  and  launch  technological 
advances.  The  decline  during  the  1970s  results  from  efforts  in  launching  heavy  payloads 
during  the  Apollo  era  and  the  race  to  the  Moon.  The  increase  in  small  satellite  launches 
during  the  1980s  and  early  1990s  stemmed  from  the  communication  industries  taking 
advantage  of  the  global  perspective  benefits  space  offered.  The  Kosmos  satellites  formed 
the  Commonwealth  of  Independent  States’  (CIS)  military  tactical  communications  con¬ 
stellation.  The  Kosmos  program  contributed  greatly  to  role  small  satellites  played  in 
the  communication  industry  during  those  two  decades.  However,  the  size  of  communi¬ 
cation  satellites  grew  as  companies  wanted  larger  systems  to  meet  increasing  customer 
needs.  The  construction  of  massive  satellites  continued  until  several  high  profile  failures 
occurred  in  the  early  1990s.  The  loss  of  the  Mars  Observer  and  the  flaw  in  the  Hubble 
Space  Telescope’s  primary  mirror  forced  government  programs  to  start  thinking  ’’faster, 
better,  cheaper”  [96] .  The  multibillion  dollar  failures  triggered  the  recent  boom  in  small 
satellites,  such  as  the  Globalstar  constellation  shown  in  Figure  1.4.  The  large  increase 
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□  Mini 
■  Micro 

□  Nano 


Figure  1.3:  Number  of  Small  Satellites  Successfully  Launched 


Table  1.1:  Mass  Classification  of  Small  Satellites 


Classification 

Wet  Mass 

Mini 

100-500  kg 

Micro 

10-100  kg 

Nano 

1-10  kg 

Pico 

<  1  kg 

5 


Figure  1.4:  The  Globalstar  Configuration 


in  small  satellite  usage  also  reflects  businesses’  and  universities’  growing  contributions 
to  the  utilization  of  space[135]. 

Today,  the  satellite  industry  is  no  longer  exclusive  to  businesses  and  government 
agencies.  Universities  are  now  using  satellites  as  educational  tools  and  platforms  for 
conducting  advancements  in  research  (see  Figure  1.5).  The  decreased  cost  and  com¬ 
plexity  of  small  satellites  when  compared  with  conventional  satellites  are  what  make 
small  satellites  more  appealing  to  academic  institutions.  Figure  1.6  shows  how  program 
costs  for  a  conventional  satellite  add  up [235].  Component  selection  and  traceability,  for¬ 
mal  design  review  costs,  infrastructure,  and  corporate  overhead  contribute  to  inflating 
the  costs  of  the  conventional  program. 

Another  benefit  to  university  programs  is  the  reduced  development  timeline  of 
small  satellites.  The  schedule  of  a  conventional  program  (3-5  years)  is  cumbersome  in 
an  academic  environment.  A  two-year  schedule  appeals  to  an  undergraduate  institution 
where  students  complete  the  bulk  of  the  advanced  courses  in  the  final  two  years  of 
the  degree.  In  addition,  the  two-year  schedule  is  beneficial  to  graduate  programs.  The 
course  work  is  completed  in  the  first  year  and  research  credits  are  earned  in  the  following 
two  or  three  years.  Figure  1.7  illustrates  a  typical  program  timeline  for  a  small  satellite 
with  a  pre-existing  infrastructure,  supporting  faculty  numbering  12,  and  approximately 
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Figure  1.5:  Distribution  of  Small  Satellite  Users[235] 


Minimum  Cost  Program: 

Price:  $2.0  M 

5x  Higher  Pointing  Acc:  $6.1  M 

Hi-Retiability  Specification:  $7.5  M 
Full  Traceability:  $12.1  M 

Formal  Periodic  Reviews:  $14.5  M 

Large  Business  Overhead: 
Conventional  Program  Price:  $17.8  M 


Figure  1.6:  Program  Costs  for  a  Typical  Conventional  Program 
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40  students  using  similar  designs  from  past  successes [235]. 


0  3  6  9  12  15  18  21  24 

Time  (months) 

Figure  1.7:  Typical  Small  Satellite  Program  Timeline 


Reduced  budget  costs  and  compact  schedules  are  not  the  only  two  aspects  of  small 
satellites  which  are  desirable  by  academic  institutions.  Universities  will  assess  mission 
failure  differently  than  how  a  commercial  operation  will.  A  business  will  experience  a 
loss  of  investment  and  time  with  each  failed  operation.  While  experiencing  the  same 
loss,  albeit  on  a  smaller  scale,  a  university  gains  the  practical  experience  of  designing, 
building,  testing,  launching  and  operating  a  satellite.  Experience  is  an  effective  edu¬ 
cational  tool  for  students.  An  orbital  failure  of  a  satellite  still  results  in  an  academic 
success  for  students.  One  of  the  mission  statements  from  the  United  States  Air  Force 
Academy’s  (USAFA)  Department  of  Astronautics,  located  in  Colorado  Springs,  CO, 
is  to  ’’teach  space  by  doing  space.”  Cadets  that  work  in  the  Space  Systems  Research 
Center  (SSRC)  are  engaged  in  each  phase  of  a  small  satellite’s  evolution.  A  corporate 
engineer  designing  a  similar  satellite  may  only  participate  in  a  single  phase  of  satellite 
development,  such  as  the  testing  or  operational  phase. 

Although  small  satellites  provide  the  aforementioned  advantages  over  their  larger 


counterparts,  several  disadvantages  exist  as  well.  The  compactness  of  a  small  satellite 
limits  the  mass  and  volume  of  usable  components.  Reduced  availability  of  space  may  ex¬ 
clude  the  use  of  more  sophisticated  components.  Compactness  may  also  greatly  increase 
the  acquisition  costs  of  miniaturized  devices  to  produce  similar  results.  In  addition,  the 
assembly  of  miniature  components  in  a  confined  space  presents  difficulties  with  arrang¬ 
ing  and  connecting  components.  Smaller  exterior  surface  area  of  the  satellite  reduces 
the  size  of  solar  arrays.  Thus,  power  supplied  to  components  is  limited.  Power  budgets, 
as  well  as  efficient  power  usage,  are  becoming  increasingly  important.  Certain  aspects  of 
operations,  such  as  payload  functionality  or  attitude  maneuver  sequences,  may  become 
impaired  if  component  power  consumption  is  greater  than  expected.  Another  strong 
disadvantage  of  a  small  satellite  is  the  difficulty  in  providing  redundant  systems  within 
the  satellite.  While  redundancy  may  be  provided  for  mission  critical  items  such  as  CPU 
memory  or  communication  systems,  it  is  often  difficult  to  incorporate  back  up  systems 
for  attitude  control  of  the  spacecraft.  If  a  component  ceases  to  perform  as  expected, 
it  may  impair  the  mission  or  cause  it  to  fail  all  together.  Therefore,  a  small  satellite 
must  rely  on  systems  that  are  not  ’power  hungry’,  do  not  occupy  too  much  space,  or 
are  fairly  simplistic  and  reliable  to  mitigate  risk. 

1.2  Gravity  Gradient  Stabilization 

The  Attitude  Determination  and  Control  Subsystem  (ADCS)  has  traditionally 
been  too  large  and  expensive  for  use  in  small  satellites.  Key  aspects  of  payload  require¬ 
ments  depend  upon  strict  pointing  requirements  and  attitude  knowledge.  The  payload 
mission  may  be  jeopardized  if  the  ADCS  subsystem  fails  in  meeting  pointing  and  atti¬ 
tude  knowledge  requirements.  To  reduce  the  probability  of  a  satellite  failing  its  mission, 
redundant  and  sophisticated  components  are  utilized  to  ensure  the  ADCS  is  not  a  sin¬ 
gle  point  of  failure  for  the  satellite.  This  poses  a  problem  for  small  satellites  requiring 
attitude  control.  In  fact,  a  common  approach  with  early  small  satellites  was  to  provide 
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Table  1.2:  Attitude  Control  Actuator  Accuracies  [69] 


Actuator 

Classification 

Accuracy 

Spin 

Passive 

±0.1° 

Magnetic  Rods 

Passive 

±5° 

Gravity  Gradient 

Passive 

±5° 

Nutation  Damper 

Passive 

±3° 

Magnetic  Torquers 

Active 

±0.1° 

Thrusters 

Active 

±0.1° 

Reaction/Momentum  Wheels 

Active 

±0.01° 

Control  Moment  Gyros 

Active 

±0.001° 

no  stabilization  at  all. 

Passive  stabilization  techniques  offered  a  way  of  providing  some  control  while 
staying  in  compliance  with  the  satellite’s  power  and  mass  requirements.  Spin,  magnetic 
rods,  and  gravity  gradient  fall  into  the  passive  category  since  no  energy  or  control 
actuators  are  required  for  operation  (although  some  stored  energy  is  required  during 
the  initial  spin- up  phase  or  deployment  of  the  boom).  These  systems  operate  in  open 
loop  (i.e.  no  control  feedback  information  is  provided  from  attitude  sensors).  Once 
the  satellite  is  in  the  desired  attitude,  it  will  exhibit  typical  pointing  performance  of 
±10°.  Active  control  systems,  (i.e.  magnetorquers,  reaction  wheels,  momentum  wheels, 
thrusters),  employ  actuators  to  generate  a  control  torque  to  the  satellite.  Active  systems 
require  attitude  sensors  to  provide  position  and  rate  information  to  close  the  loop  of 
the  control  system.  By  using  feedback  information  to  fine  tune  the  control  system, 
performance  values  of  active  control  will  be  at  least  an  order  of  magnitude  smaller  than 
passive  systems.  Table  1.2  lists  various  attitude  control  actuators,  their  classification  as 
passive  or  active,  and  typical  accuracy  values. 

While  passive  systems  were  used  in  the  early  prototypes  of  small  satellites  (be¬ 
tween  1957  and  1980),  active  systems  are  seeing  more  use  in  supplementing  passive 
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Figure  1.8:  Structural  Drawing  of  UoSAT-1 


components.  When  a  gravity  gradient  boom  is  carried,  a  pitch  wheel  can  be  used  in 
momentum  bias  mode  or  a  yaw  wheel  in  zero-momentum  bias.  The  most  common  con¬ 
trol  system  for  current  small  satellites  is  a  combination  of  gravity  gradient  boom  and 
magnetorquers.  This  control  scheme  allows  the  satellite  to  have  yaw  pointing  control 
while  remaining  nadir  pointing.  The  first  small  satellite  to  use  this  control  configuration 
was  UoSAT-1  (OSCAR-9),  an  amateur  satellite  built  and  operated  by  the  University  of 
Surrey  in  the  early  1980s  (see  Figure  1.8)  [218]. 

Gravity  gradient  stabilization  is  not  a  new  concept.  D’Alembert  and  Euler’s  ce¬ 
lestial  mechanics  work  in  1749  first  discussed  gravitational  gradient  effects  of  an  axially 
symmetric  ellipsoid  in  an  inverse  square  field.  30  years  later  Lagrange  used  this  infor¬ 
mation  to  explain  the  librations  of  the  Moon[159].  Uneven  mass  distribution  within  the 
Moon’s  crust  results  in  a  gravity  gradient  stabilized  attitude  with  respect  to  the  Earth. 
Astronomers  have  observed  this  gravity  gradient  stabilized  attitude  for  many  years  and 
have  termed  the  unseen  50%  as  the  ’dark  side’  of  the  Moon.  Actually,  due  to  the  lon¬ 
gitudinal,  latitudinal,  and  diurnal  components  of  the  Moon’s  apparent  libration,  over  a 
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Figure  1.9:  Gravity  Gradient  Stabilization  with  a  Momentum  Wheel 


period  of  time  an  Earth  based  observer  is  able  to  view  59%  of  the  Moon’s  surface[176]. 

Gravity  gradient  stabilization  of  a  satellite  occurs  when  the  minimum-inertia  axis 
locks  in  a  vertical  orientation  with  the  orbiting  body.  This  orientation  is  conventionally 
accomplished  through  the  use  of  a  deployable  boom  with  a  tip  mass  at  the  end.  Once 
the  boom  is  deployed,  the  tip  mass  will  increase  the  moment  of  inertia  in  the  directions 
transverse  to  the  boom[240].  This  boom  configuration  allows  passive  control  in  both 
the  roll  and  pitch  directions;  therefore  a  communication  or  remote  sensing  payload 
may  remain  nadir  pointing.  This  boom  configuration  does  not,  however,  provide  yaw 
control  since  the  satellite  is  free  to  rotate  about  the  vertical  axis.  Early  research  efforts 
illustrated  how  a  stable  orientation  is  created  when  a  momentum  wheel  with  its  angular 
momentum  is  aligned  along  the  positive  orbit  normal  (shown  in  Figure  1.9)  [223].  A 
detailed  discussion  of  the  theoretical  development  for  gravity  gradient  stabilization  is 
provided  in  Appendix  A  of  this  document. 

Fairing  limitations  of  launch  vehicles  requires  the  containment  of  boom  systems 
prior  to  launch.  The  basic  idea  of  deploying  gravity  gradient  booms  while  in  orbit 
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was  initially  developed  and  demonstrated  in  space  as  early  as  the  1960’s[160].  The 
deployment  process  involves  storing  strain  energy  within  a  compact  structure  prior  to 
fitting  the  satellite  to  the  launch  vehicle.  Once  the  satellite  is  on-orbit,  an  initiation 
command  is  given  by  either  ground  based  operators  or  an  autonomous  control  system 
contained  within  the  satellite’s  data  handling  system.  This  initiation  command  signals 
the  gravity  gradient  component  to  release  the  stored  strain  energy  which  deploys  the 
boom  element  to  obtain  its  final,  desired  shape  or  configuration.  Conventional  boom 
systems  use  telescopic  sections  (shown  in  Figure  1.10)  or  furlable  overlapped  tubing, 
often  referred  to  as  Storable  Tubular  Extendible  Member  (STEM),  as  shown  in  Fig¬ 
ure  1.11.  The  telescopic  boom  built  by  SSTL,  known  as  STACER,  has  a  heritage  of 
over  25  years.  During  this  time,  spacecraft  and  sounding  rockets  have  utilized  more  than 
600  STACER  units[135].  Small  satellite  designers  find  STACER’s  compact  dimensions 
while  stowed  (102  x  115  x  264  mm)  and  lightweight  design  (2.2  kg  without  tip  mass) 
very  appealing.  STACER’s  telescopic  section  is  spring  loaded.  Deployment  is  initiated 
with  dual  redundant  pyrotechnic  bolt  cutters.  The  ” pyro-cutters” ,  classified  as  ’’Class 
C”  explosives,  cut  through  a  shear  bolt  to  release  the  stored  strain  energy. 

STEM  booms  use  heat  treated  steel  to  deploy  a  rigidizable  antenna.  The  boom 
element  is  fabricated  to  a  specified  length  and  rolled  like  a  sleeping  bag.  For  longer 
lengths,  it  is  divided  into  segments  that  are  joined  by  a  thin  lap  joint.  More  recently, 
AEC-Able  Engineering  developed  the  CoilAble  boom  (shown  in  Figure  1.12).  The 
CoilAble  boom  is  fabricated  as  a  helix  and  stored  in  a  collapsed  configuration  inside  the 
deployment  canister [116].  However,  STEM’s  tend  to  be  heavy  due  to  the  use  of  either 
beryllium  copper  or  stainless  steel. 

Although  these  conventional  booms  do  have  significant  flight  heritage,  they  pos¬ 
sess  several  unfavorable  characteristics.  The  high  amount  of  stored  energy  in  springs  and 
the  use  of  pyrotechnic  bolt  cutters  generate  survivability  concerns  for  a  launch  vehicle. 
Elaborate  inhibit  measures  need  to  be  applied  during  launch  to  ensure  an  early  release 
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Figure  1.10:  Illustration  of  the  STACER  Telescopic  Section  Boom 


EXTENSION  MECHANISM 


Figure  1.11:  Example  of  a  STEM  Boom  [67] 


Figure  1.12:  AEC-Able  Engineering’s  CoilAble  Boom  [49] 
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does  not  occur.  An  early  release  could  prove  fatal  not  only  to  the  satellite  and  launch 
vehicle,  but  to  other  satellites  sharing  a  ride  on  the  same  launch  vehicle  as  well.  Also, 
human  life  may  be  at  risk  if  the  Space  Shuttle  is  used  to  launch  a  satellite  that  utilizes 
conventional  type  boom  deployment  systems.  The  STACER  boom  design  suffers  from 
high  packaged  strain  energy.  This  packaged  strain  energy  limits  the  size  of  the  boom  to  a 
maximum  diameter  of  2”  and  produces  low  deployed  stiffness  and  strength  [36].  In  addi¬ 
tion,  pyrotechnic  release  devices  produce  high  shock  and  contamination,  are  hazardous 
and  costly  to  handle,  and  are  not  re-setable. 

One  alternative  to  using  a  traditional  boom  is  to  inflate  a  structure  while  on  orbit. 
Deployment  tests  of  inflatable  booms  have  been  extensively  conducted  during  the  last 
fifteen  years.  Haug  and  colleagues [85]  utilized  finite  element  methods  to  simulate  the 
deployment  process  of  an  inflatable  antenna  in  1991.  Six  years  later,  Tsoi’s  thesis[226] 
used  these  equations  to  simulate  the  deployment  of  inflatable  booms  which  were  both 
folded  and  rolled  up.  Clem  and  his  associates [46]  conducted  deployment  tests  in  2001 
while  a  year  later  Campbell  and  her  coworkers [35]  experimentally  investigated  gravita¬ 
tional  effects  during  the  deployment  of  inflatable  tubes.  Also  during  2002,  Miyazaki  and 
Uchiki[153]  and  Wang  and  Johnson[232]  tested  and  analyzed  the  deployment  process  of 
inflatable  structures.  All  of  these  studies  focused  on  pressure  stabilized  inflatable  booms 
where  constant  pressure  is  required  to  maintain  the  rigidity  of  the  structures. 

While  the  mentioned  studies  on  inflatable  structures  were  thorough,  they  only 
focused  on  the  deployment  dynamics  with  minimal  consideration  given  to  the  control 
authority  impact  on  the  spacecraft.  In  addition,  the  use  of  inflatable  structures  were 
never  shown  to  be  a  feasible  system  for  small  satellites.  The  complexity  of  the  plumbing 
system  for  the  pressurized  air  as  well  as  the  added  mass  and  volume  requirements 
for  storage  tanks  would  quickly  exceed  the  constraints  placed  upon  small  satellites. 
Therefore,  the  support  structure  needed  for  inflatable  structures  would  be  scaled  down 
to  fit  within  the  confines  of  a  small  satellite  and  greatly  reduce  the  length  of  such  a 
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structure. 

1.3  Elastic  Memory  Composites 

A  recent  innovation  which  has  the  potential  to  revolutionize  the  design  of  de¬ 
ployable  booms  is  the  development  of  shape  memory  materials  and  their  introduction 
into  the  design  of  deployment  mechanisms.  Shape  memory  ’’mechanisms”  can  elimi¬ 
nate  the  need  for  traditional  highly  complex  mechanical  deployment  devices,  massive 
launch  canisters,  and  independent  deployment  control  systems.  In  addition,  these  shape 
memory  mechanisms  can  lead  to  dramatically  simpler  boom  designs  that  include  fewer 
’’parasitic”  (i.e.  non-structural)  parts  and  are  therefore  much  lighter  in  weight. 

A  shape-memory  alloy  (SMA)  is  a  mixture  of  two  or  more  metals  that  has  the 
special  property  of  being  able  to  ’’memorize”  a  certain  shape,  and  return  to  that  shape 
even  after  being  deformed.  Usually  the  return  to  a  memorized  shape  is  triggered  by 
heating  the  material.  The  change  can  happen  very  fast,  often  within  seconds.  The  key 
to  a  SMA’s  ability  to  shape-change  is  that  its  structure  differs  depending  on  the  temper¬ 
ature.  At  high  temperatures,  the  atoms  in  an  SMA  possess  a  very  stiff,  rigid  structure, 
called  the  ’’austenitic”  structure  (named  after  British  metallurgist  Sir  William  Chandler 
Roberts- Austen).  The  shape  of  an  SMA  is  linked  to  its  austenitic  structure [225].  Any 
change  in  the  shape  of  the  metal  while  it  is  in  the  austenitic  phase  causes  the  structure 
to  change,  and  vice  versa. 

As  the  metal  cools  and  reaches  a  critical  temperature  range,  the  atoms  begin  to 
realign  themselves  into  a  different  structure,  called  the  ’’martensitic”  structure  (named 
after  German  metallurgist  Adolf  Martens).  This  structure  is  also  linked  to  the  austenitic 
structure,  but  is  flexible  and  allows  the  metal  to  be  visibly  bent,  stretched,  and  manipu¬ 
lated  without  changing  the  underlying  atomic  structure.  When  the  metal  is  heated  again 
to  its  critical  temperature  range,  the  metal  transforms  back  into  the  rigid  austenitic 
structure  which,  being  linked  to  the  shape  of  the  metal,  causes  it  to  regain  its  original, 
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’’memorized”  or  ’’programmed”  shape[73].  If  the  metal  is  cooled  once  again  to  the  flex¬ 
ible  martensitic  phase,  it  will  retain  the  memorized  shape  until  otherwise  changed  by 
an  outside  influence. 

The  austenitic  and  martensitic  structures  that  cause  shape  memory  also  possess 
several  other  exceptional  properties:  superelasticity,  damping,  and  high,  controllable 
recovery  force.  The  superelasticity  of  an  SMA  enables  the  SMA  to  stretch  more  than 
other  metals  and  then  spring  back  to  original  size.  SMAs  can  be  stretched  up  to  8%  more 
than  their  original  length  without  permanently  stretching  or  damaging  the  material 
(other  metals  stretch  less  than  1%). 

Damping,  another  property  of  shape  memory  composites,  is  the  ability  to  stop 
oscillations  or  vibrations  quickly.  For  example,  when  a  rubber  ball  is  dropped  onto  a 
firm  surface,  it  will  bounce  several  times,  but  each  successive  bounce  will  get  lower  and 
lower,  until  eventually  the  ball  stops  bouncing.  On  a  concrete  surface  the  ball  may 
bounce  many  times,  whereas  the  ball  might  only  bounce  once  or  twice  on  a  carpeted 
surface.  The  carpeted  surface,  then,  would  have  high  damping  (it  stops  the  ball  quickly), 
and  the  concrete  surface  has  lower  damping.  The  same  idea  applies  to  damping  in  a 
mechanical  system.  An  automotive’s  engine  produces  numerous  vibrational  frequencies. 
The  amplitude  of  these  vibrations  experienced  by  the  occupants  depends  upon  the 
damping  capability  of  the  car’s  materials  and  structure.  Since  SMAs  have  high  damping, 
a  person  sitting  in  an  operating  SMA  car  might  not  feel  the  engine  vibrating. 

A  high  recovery  force  from  an  SMA  can  be  compared  to  the  stretching  or  com¬ 
pression  of  a  spring,  such  as  one  found  in  a  mattress  or  a  trampoline.  When  a  person 
jumps  onto  a  trampoline,  the  springs  are  compressed.  When  the  springs  return  to  their 
normal  state,  the  restoring  force  sends  the  person  flying  into  the  air.  Changing  the 
shape  of  an  SMA  in  its  martensitic  phase  is  like  compressing  the  springs  on  a  tram¬ 
poline.  The  stored  strain  energy  is  released  when  the  SMA  returns  to  its  memorized 
austenitic  shape. 
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Other  shape-memory  materials  also  have  these  properties  (damping,  superelas¬ 
ticity,  etc).  Shape-memory  polymers  (SMPs)  are  a  type  of  shape- memory  material  that 
is  also  triggered  by  temperature  change,  but  is  composed  of  plastic,  instead  of  metal. 
The  most  familiar  example  of  an  SMP  is  shrink-wrap [11].  SMPs  possess  greater  levels 
of  superelasticity  than  SMAs.  SMPs  can  stretch  up  to  400%  more  than  their  original 
length.  However,  SMPs  produce  a  lower  recovery  force.  This  means  SMPs  are  more  fa¬ 
vorable  for  purposes  where  major  shape  changing  needs  to  occur  while  generating  little 
resistance  to  the  change.  An  example  usage  of  SMPs  which  is  currently  being  mar¬ 
keted  is  ’’smart”  cold- weather  clothing.  The  clothing  is  impermeable  to  wind  and  rain 
at  lower  temperatures  but  porous  at  higher  temperatures  for  ’’breathability.”  Another 
difference  between  SMPs  and  SMAs  is  that  SMPs  can  be  manufactured  from  biodegrad¬ 
able  materials.  This  is  advantageous  for  surgical  purposes.  It  is  not  necessary  to  remove 
temporary  implants  or  sutures  created  from  SMP  materials  since  they  will  intentionally 
disintegrate  over  time[122]. 

The  first  use  of  SMAs  occurred  in  various  fields  of  engineering,  especially  in 
the  military.  The  first  SMAs  were  developed  by  the  Navy.  A  nickel-titanium  alloy 
was  developed  at  the  US  Naval  Ordnance  Laboratory  in  1961.  The  shape-memory 
properties  were  discovered  accidentally.  The  Navy  named  the  new  substance  ’’Nitinol”— 
Ni  for  nickel,  Ti  for  titanium,  and  NOL  for  Naval  Ordnance  Lab[105].  Nitinol  has  been 
the  most  prominently  used  SMA  in  engineering  applications. 

The  first  application  of  SMAs  was  for  pipe  joining.  For  this  purpose,  the  SMA 
is  programmed  into  the  shape  of  a  short  tube  slightly  smaller  than  the  two  pipes  being 
joined.  Then,  in  its  malleable  phase,  the  metal  is  stretched  into  a  tube  slightly  larger 
than  the  two  pipes.  The  SMA  is  set  with  the  two  pipe  ends  inside,  and  then  warmed. 
As  warming  occurs,  the  pipes  are  squeezed  together  and  secured.  This  method  has  been 
used  quite  successfully  to  join  various  types  of  pipes,  from  hydraulic  pipes  in  F-14  fighter 
planes[73]  and  naval  ships,  to  transport  pipes  in  the  chemical  and  petroleum  industries. 
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Pipes  of  up  to  80  inches  in  diameter  have  been  joined  utilizing  SMAs. 

Due  to  the  controlled  force  exerted  when  an  SMA  regains  its  programmed  shape, 
SMAs  can  be  used  in  robotics  to  provide  smooth  motions.  Robotics  engineers  in  Japan 
have  developed  ’’devices  to  grasp  delicate  paper  cups  filled  with  water”  [178].  SMAs  have 
also  been  used  as  replacements  for  explosives  in  building  demolitions.  When  a  large 
force  stored  in  an  SMA  is  recovered  in  fractions  of  a  second,  the  yield  is  comparable 
to  an  explosion.  These  two  examples  of  the  practical  uses  of  SMAs  demonstrate  these 
materials’  astonishing  range  of  capabilities. 

SMAs  are  frequently  utilized  in  the  aerospace  industry.  The  high  damping  capa¬ 
bility  of  SMAs  provides  a  spacecraft  insolation  from  vibration,  a  major  concern  during 
launch.  In  addition,  SMAs  can  be  utilized  for  deploying  instruments  and  payloads 
once  the  spacecraft  reaches  orbit.  Many  tasks,  from  the  simple  uncovering  of  a  camera 
lens,  to  the  complete  deployment  of  a  solar  panel,  can  be  activated  by  an  SMA.  These 
actions  can  be  programmed  to  trigger  automatically  once  the  alloy  reaches  a  prede¬ 
termined  low  temperature  as  it  cools  off  in  the  space  environment.  SMAs  can  also  be 
activated  through  heating  from  an  electrical  current.  The  heating  process  can  be  easily 
commanded  from  ground  control  at  the  appropriate  time. 

Elastic  Memory  Composite  (EMC)  materials  are  a  relatively  new  addition  to  the 
family  of  SMAs  and  have  been  under  development  by  Composite  Technology  Devel¬ 
opment,  Inc.  (CTD)  since  1999  [227].  CTD  has  been  able  to  dramatically  improve 
both  the  stiffness  and  the  recovery  force  of  a  shape  memory  polymer  by  incorporating 
it  into  a  fiber-reinforced  composite[125].  This  results  in  substantially  lower  densities 
and  higher  elastic  strain  capacities.  A  specific  thermo-mechanical  cycle  is  used  to  store 
and  release  the  strain  within  the  material.  The  strains  are  induced  by  elevating  the 
temperature  of  the  material  above  the  glass  transition  temperature  (Tg)  of  the  resin 
and  then  applying  a  mechanical  load  to  deform  the  material.  Venturing  above  this  crit¬ 
ical  temperature  results  in  a  ’rubbery’  composite  structure  which  allows  high  levels  of 
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strain  to  be  achieved  without  damaging  the  composite  fibers  or  the  resin  system.  While 
maintaining  the  mechanical  load,  the  material  is  cooled  below  Tg.  This  ’freezes’  the 
strain  within  the  material.  Now,  the  mechanical  load  can  be  removed  and  the  EMC 
will  remain  deformed  until  it  is  once  again  heated  above  Tg.  The  heated  material  will 
release  the  strain  energy  stored  within  and  will  return  to  its  original  shape [71]. 

A  study  was  conducted  in  2002  by  members  of  CTD  and  ABLE  Engineering[115]. 
The  purpose  of  the  study  was  to  determine  the  impact  of  using  EMC  longerons  in  the 
CoilAbleTM  boom  design  for  NASA’s  Space  Solar  Panel  System.  The  study  determined 
that  the  longeron  mass  of  the  boom  could  be  reduced  up  to  a  factor  of  ten  while  easily 
achieving  effective  strains  of  2%.  Repeated  stowage  and  deployment  cycles  showed  no 
substantial  mechanical  degradation. 

1.4  Flexible  Spacecraft  Control 

With  the  conclusion  of  the  Apollo  space  program  in  the  1970s,  NASA  researchers 
turned  their  attention  to  the  design  of  flexible  space  structures.  This  generated  an 
increasing  interest  in  the  dynamics  and  control  of  such  structures  within  various  re¬ 
search  communities.  Several  survey  papers  exist  which  provide  a  listing  of  published 
efforts  in  the  area  of  controlling  flexible  structures.  Robinson’s  survey[177]  lists  a  small 
number  of  papers  on  structural  control,  and  in  particular  on  attitude  control  of  flex¬ 
ible  space  structures.  Croopnick  et  al.  [50]  present  a  literature  survey  in  the  areas  of 
attitude  control,  vibration  control  and  shape  control  as  they  apply  to  space  structures. 
Meirovitch[149]  assesses  various  methods  for  the  active  control  of  space  structures  with 
a  view  to  the  problems  of  high  dimensionality  and  modeling.  Balas[15]  presents  a  math¬ 
ematical  framework  for  the  discussion  of  large  space  structure  (LSS)  control  theory  and 
provides  a  look  at  trends  in  LSS  control  theory  in  the  early  1980s.  A  comprehensive 
survey  of  problems  in  dynamic  modeling  and  control  of  space  structures  is  compiled 
by  Nurre  et  al.  [161] .  One  important  thing  to  note  from  the  surveys  listed  above  is  the 
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focus  on  the  control  of  large  flexible  space  structures;  a  focus  which  has  predominately 
driven  the  field  of  study. 

Flexible  structures  on  small  satellites  are  becoming  increasingly  popular.  Actively 
controlling  flexible  structures  is  an  important  area  of  research  [11 7].  The  increased  usage 
of  flexible  structures  is  related  to  the  growing  number  of  small  satellites  being  launched, 
the  use  of  lightweight  materials,  and  the  desire  to  meet  tighter  packaging  requirements 
inside  a  launch  vehicle’s  fairing[133].  Tethering  spacecraft,  deployable  solar  arrays, 
communication  antennas,  and  extending  scientific  payloads  away  from  the  spacecraft 
bus  are  only  a  few  examples  where  pointing  accuracies,  attitude  information,  and  active 
control  are  required.  The  excitation  of  these  appendages  is  highly  probable  while  the 
spacecraft  performs  attitude  control  maneuvers.  Reorientation  of  the  spacecraft  during 
these  maneuvers  is  a  critical  part  of  the  ADCS  process  so  that  the  induced  vibrations 
are  kept  to  minimum  levels  [77]. 

Maghami,  Sparks,  and  Lim  (1998)  identified  four  dynamic  characteristics  of  flex¬ 
ible  structures  which  complicate  control  system  design[142]. 

•  large  number  of  structural  modes  in  controller  bandwidth 

•  low  and  closely  spaced  modal  frequencies 

•  very  small  inherent  damping 

•  insufficient  knowledge  of  parameters 

Even  if  the  underlying  physical  model  of  the  spacecraft  could  be  accurately  modeled  at 
one  point  in  time,  parameter  variations  during  system  operation  will  eventually  result 
in  an  inaccurate  model[236].  The  distributed  parameter  models  often  used  to  describe 
flexible  structures  are  essentially  infinitely  dimensional  and  need  to  be  truncated.  This 
leads  to  complications  of  destabilizing  one  or  more  of  the  poorly  damped  modes  since 
only  some  of  the  lower  frequency  modes  are  approximated  and  kept  [130]. 
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Research  in  the  area  of  solid-state  physics  originally  revealed  how  modal  distortion 
may  impede  the  propagation  of  vibrational  energy  [9].  Anderson’s  efforts  led  to  improved 
predictions  on  the  effects  defects  have  in  lattice  vibrations.  Almost  five  decades  later, 
the  theories  derived  from  solid-state  physics  are  being  applied  to  flexible  spacecraft 
structures [246].  When  modal  distortions  occur,  they  tend  to  be  localized  in  the  region 
in  which  they  were  created.  Most  spacecraft  flexible  appendages  are  unactuated  or 
unsensed.  Slight  structural  asymmetries  may  complicate  control  design  as  the  potential 
for  large  modal  errors  compared  to  modeled  results  will  push  for  an  overly  conservative 
control  design[150]. 

Another  technique  sensitive  to  modeling  errors  was  the  wave  cancelation  technique  [203]. 
This  technique  drove  a  second-order  system  to  its  final  position  in  finite  time.  In  the 
early  1990s,  researchers  furthered  this  technique  using  a  pulsed  sequence  expansion  to 
desensitize  their  models[197].  Input  shaping  relies  on  the  fact  that  a  spectrum  of  a 
convolution  of  two  signals  is  a  product  of  each  signal’s  spectra.  Thus,  a  zero  excitation 
frequency  results  when  one  or  the  other  spectra  is  zero.  By  designing  a  sequence  of 
input  pulses,  Singer  and  Seering  (1990)  forced  the  magnitude  of  the  residual  energy  to 
zero.  A  time  delay  filter  designed  to  cancel  the  poles  of  the  system  produced  similar 
results  while  remaining  insensitive  to  errors  in  modeled  damping  and  frequency  [199]. 
Multiple  zeros  of  the  time  delay  filter  could  also  be  placed  at  the  estimated  locations  of 
the  system  poles  to  produce  robust  time-optimal  control[131]. 

The  general  problem  with  frequency-based  input  shaping  is  they  are  designed 
for  linear  systems.  A  small  satellite  with  extended  flexible  appendages  is  subjected  to 
many  disturbance  torques,  several  of  which  are  unmodeled.  If  the  input  shaping  method 
accounts  for  nonlinearities  by  making  robustness  assumptions,  the  spacecraft  bus  may 
experience  significant  levels  of  residual  vibration[107].  Time-domain  representations 
and  optimization  of  control  inputs  are  usually  easier  to  generalize.  Non  frequency-based 
input  shaping  techniques  use  either  some  type  of  inverse  dynamics  computation  with  a 
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known  endpoint  trajectory  or  an  optimization  technique  to  find  a  feedforward  function 
to  meet  final  conditions[170].  Some  additional  types  of  input  shaping  techniques  either 
increase  the  range  of  uncertain  parameters  in  regions  where  residual  vibration  is  below 
a  specified  level [200]  or  include  probability  distributions  of  the  uncertain  parameters  to 
weight  the  nominal  value  of  the  parameter  the  most  [166]. 

Input  shapes  are  designed  to  result  in  an  exact  time-optimal  control  algorithm. 
To  do  this,  a  ’’bang-bang”  method  of  actuation  is  used.  This  method  is  an  initial  firing  of 
the  control  actuators  and  is  used  to  begin  a  maneuver;  then  the  spacecraft  is  placed  into 
a  ’coast’  mode.  When  the  attitude  sensors  detect  the  system  approaching  the  opposite 
threshold,  the  actuators  fire  in  the  opposite  direction  to  stop  the  momentum  from  the 
initial  maneuver.  This  attitude  maneuver  process  continues,  while  the  actuator  inputs 
decrease  in  magnitude,  until  a  narrow  band  of  motion  is  achieved  about  the  desired  final 
position.  The  ADCS  then  monitors  the  satellite’s  motion  until  once  again  the  sensors 
determine  an  attitude  threshold  is  about  to  be  breached. 

The  ’’bang-bang”  method  attempts  to  maximize  the  number  of  zero  excitation 
frequencies  in  the  input  spectra.  However,  higher-order  frequencies  of  the  system  can 
easily  be  excited  with  this  method.  A  smoother  control  input  is  achieved  by  using  an 
approximation  function  to  eliminate  the  sudden  change  of  control  magnitudes  [102].  The 
smoother  control  input  is  less  likely  to  excite  the  structural  modes  and  is  calculated  by 
solving  an  optimal  control  problem  with  the  objective  of  minimizing  the  maneuver  time. 
The  solution  utilizes  the  state  equations  describing  the  rigid-body  mode,  boundary  con¬ 
ditions,  and  an  additional  constraint  which  limits  the  derivative  of  the  control  input[7]. 
Numerical  solutions  of  Albassam’s  (2002)  technique  show  the  residual  energy  of  the 
flexible  appendages  is  greatly  reduced  while  slightly  increasing  the  maneuver  time. 

Accounting  for  modeling  errors  and  uncertainties  in  control  methods  have  been 
studied  by  several  researchers  [98]  [99]  [108]  [109]  [146]  [187]  [188]  [189]  [190]  [201]  [202]  [243]  [244] 
However,  one  paper  is  of  particular  interest.  Doyle  and  Stein [58]  investigated  the  is- 
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sue  of  feedback  control  design  in  the  face  of  uncertainties  and  generalized  single-input, 
single-output  statements  and  constraints  of  the  design  problem  to  multi-input,  multi¬ 
output  cases.  They  developed  the  procedure  known  as  linear-quadratic-Gaussian/loop 
transfer  recovery  (LQG/LTR).  Stein  and  Athans[213]  provide  a  tutorial  overview  of  the 
LQG/LTR  procedure  for  linear  multivariable  feedback  control  design.  This  method  is 
applied  to  robust  controller  synthesis  for  a  large  flexible  space  antenna  by  Joshi[98]  and 
Sundararajan[216]  while  a  modified  version  is  presented  by  Blelloch  and  Mingori[21]. 

The  linear-quadratic-regulator  (LQR)  is  a  special  optimal  controller  whose  cost 
function,  measure  of  performance,  is  a  quadratic  function  of  states  and  controls.  Two 
desirable  properties  of  the  LQR  are  good  stability  margins  and  sensitivity  properties. 
One  limitation  is  that  it  is  a  full  state  feedback  type  of  controller.  In  several  practical 
applications,  access  to  all  of  the  states  is  difficult  to  achieve  and  state  estimation  will  be 
required.  The  LQG  problem  combines  the  LQR  controller  with  an  estimation  filter  (i.e. 
Kalman  filter).  However,  the  LQG  controller  will  often  have  lower  stability  margins  and 
lower  gain  crossover  frequency  than  the  LQR  controller.  The  LQG  will  pass  more  noise 
into  the  system  and  have  a  slower  response  when  compared  to  the  LQR. 

The  main  problem  with  the  LQG  solution  is  its  lack  of  robustness  which  has 
resulted  in  a  failure  to  work  effectively  in  real  environments[212].  As  more  realism  is 
added  to  the  plant  of  the  system,  the  LQG  became  unstable  in  the  presence  of  model 
uncertainties.  The  loop  transfer  recovery  technique  (LQG/LTR)  addresses  some  of 
the  shortcomings  of  the  LQG  approach.  The  process  begins  with  selecting  the  LQR 
parameters  until  a  desired  open-loop  transfer  function  is  obtained.  The  filter  design 
parameters  are  iterated  until  the  desired  loop  transfer  function  shape  has  been  obtained. 
This  allows  the  LQG  technique  to  become  a  flexible  frequency  domain  design  technique. 
It  is  still  based  on  state-space  design  techniques,  but  is  more  of  a  classical  approach  to 
controller  design. 

Usage  of  LQG/LTR  techniques  is  applicable  to  current  research  efforts.  In  par- 
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ticular,  Mackison[139][140][141]  has  used  these  techniques  to  develop  controllers  which 
retain  the  characteristics  of  full  state  feedback  without  requiring  all  of  the  states  to 
be  measureable.  This  work  has  demonstrated  full  state  feedback,  LQG/LTR  feedback, 
and  reduced  order  compensators  derived  from  the  LQG/LTR  controller  all  produced 
identical  dynamic  responses  for  3-axis  attitude  control  problems.  Smaller  initial  control 
torques  were  required  for  the  LQG/LTR  and  reduced  order  controllers  when  compared 
to  that  required  by  full  state  feedback  controllers.  This  result  occurs  because  a  low  order 
transfer  function  for  the  controller  is  found  through  pole-zero  cancellation  techniques. 

1.5  Problem  Formulation 

The  support  structure  size  of  satellites  are  limited  to  the  dimensions  of  the  launch 
vehicle  payload  fairing  being  used  to  deliver  the  satellite  to  orbit.  An  exception  to  this 
occurs  when  the  spacecraft  is  constructed  in  orbit.  The  International  Space  Station 
is  one  example  of  an  operational  space  structure  which  was  assembled  after  launch. 
However,  the  dimensions  of  the  space  station  sections  were  limited  by  the  size  of  the 
Space  Shuttle  delivery  bay.  A  popular  method  used  to  alter  the  support  structure  size 
of  a  satellite  is  to  package  components  prior  to  launch  and  then  deploy  these  structures 
once  the  spacecraft  is  in  orbit. 

Deployable  spacecraft  structures  are  often  used  for  communication  antennas,  cre¬ 
ating  larger  surface  areas  for  solar  arrays,  moving  scientific  payloads  away  from  the 
satellite,  and  providing  passive  gravity  gradient  stabilization.  These  structures  extend 
out  from  their  stowed  configuration  using  various  forms  of  deployment  mechanisms  and 
energy  sources.  Reliability,  low  mass,  packaged  volume,  and  energy  consumption  are 
design  concerns  for  deployable  appendages.  Since  it  is  ideal  to  make  these  structures 
mechanically  simple  and  light  weight,  they  are  susceptible  to  vibrations;  either  induced 
by  attitude  maneuvers  of  the  spacecraft  or  external  disturbance  torques.  These  vibra¬ 
tions  may  prove  detrimental  to  mission  completion  and  become  more  of  a  concern  as 
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the  inertial  properties  of  the  satellite  reduces. 

The  review  of  current  literature  in  the  areas  of  deployable  structures  and  small 
satellites  identifies  a  need  to  further  relevant  research.  High  levels  of  stored  energy 
in  mechanical  deployment  systems,  as  well  as  the  use  of  pyrotechnic  cutting  bolts, 
require  additional  levels  of  inhibit  measures  to  prevent  premature  release  during  launch 
environments.  High  packaged  strain  energy  of  traditional  booms  will  limit  the  diameter 
of  the  boom  while  also  producing  low  deployed  stiffness  and  strength.  In  addition,  the 
hazardous  pyrotechnic  release  devices  will  produce  high  shock  and  contamination. 

Mass  limitations  of  small  satellites  may  make  redundant  safety  measures  difficult 
to  use.  Reduced  usable  mass  will  also  impact  the  sizing  options  of  deployable  structures 
available  for  small  satellites.  As  the  length  of  the  appendage  increases,  so  will  its  mass. 
This  would  lead  to  ruling  out  the  use  of  small  satellites  to  accomplish  certain  missions 
which  would  require  longer  deployed  structures.  Some  examples  of  missions  using  long 
booms  are  interferometry,  generating  power  along  the  boom  length,  determining  plasma 
variations  local  to  the  satellite,  conducting  wake  studies  to  predict  spacecraft  charging, 
and  there  are  probably  many  more  that  we  have  not  yet  considered. 

Shape  memory  materials  show  promise  in  minimizing  the  shortcomings  in  tradi¬ 
tional  booms.  Although  shape  memory  technology  has  been  studied  over  the  last  50 
years  (e.g.  US  Naval  Ordnance  Laboratory  research  in  1961),  Elastic  Memory  Compos¬ 
ites  (EMC)  are  an  emerging  field  of  study.  These  composites  will  result  in  substantially 
lower  densities  and  higher  elastic  strain  capacities.  By  constructing  a  deployable  struc¬ 
ture  from  EMC  materials,  it  is  possible  to  eliminate  the  need  for  traditional  highly 
complex  mechanical  deployment  devices.  These  structures  contain  fewer  non-structural 
components  and  are  lighter  in  weight.  A  recent  study  completed  by  Composite  Technol¬ 
ogy  Development,  Inc  (CTD)  determined  longeron  mass  of  a  deployed  structure  could 
be  reduced  by  an  order  of  magnitude  while  achieving  effective  strains  of  2%  [11 5]. 

Constructing  flexible  structures  from  ’soft’  materials  will  introduce  more  flexible 
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modes  into  the  system  dynamic  equations  of  motion.  Control  designers  will  need  to 
concern  themselves  not  only  with  the  additional  flexible  modes,  but  also  the  lower 
frequencies  at  which  they  occur.  Research  over  the  last  three  decades  mainly  focused 
on  attitude  control  of  large  space  structures  with  flexible  appendages  because  a  limited 
number  of  bending  modes  would  form  close  to  the  control  bandwidth.  The  large  inertia 
values  of  these  spacecraft  helped  to  keep  the  control  bandwidth  low  enough  in  frequency 
so  that  typically  only  the  first  bending  mode  was  included  in  the  analysis. 

Attitude  control  authority  of  small  satellites  will  be  limited  more  from  both  in¬ 
creased  resonant  modes  at  lower  frequencies  as  well  as  lower  inertia  values  of  the  central 
body.  Because  of  the  added  risk,  small  satellite  engineers  are  more  likely  to  select  flex¬ 
ible  appendages  which  have  extensive  flight  heritage.  The  traditional  response  is  to 
design  around  the  flexible  nature  of  the  system  by  making  the  appendages  stiffer,  the 
spacecraft  inertia  larger,  or  lowering  the  attitude  controller’s  closed  loop  bandwidth. 
Engineers  from  Surrey,  a  leading  designer  of  small  satellites,  use  only  traditional  booms 
and  would  not  entertain  the  thought  of  appendages  made  from  ’soft’  materials.  This 
directly  results  from  the  uncertainties  involved  and  a  perceived  view  of  the  high  risks 
surrounding  non-traditional  booms.  A  concern  arises  when  mission  requirements  call 
for  the  use  of  a  small  satellite  as  well  as  a  flexible  non-traditional  structure. 

To  address  this  concern,  this  research  will  continue  following  the  path  defined  by 
Mackison[139][140][141].  While  the  LQG/LTR  technique  was  shown  effective  in  gener¬ 
ating  reduced  order  controllers  for  3-axis  attitude  control  of  spacecraft,  Mackison  only 
modeled  the  rigid  body  dynamics  of  the  spacecraft.  This  research  will  examine  the 
unique  characteristics  of  controllers  for  a  satellite  with  a  flexible  appendage  constructed 
from  EMC  materials.  In  particular,  this  research  will  address  the  following  basic  ques¬ 
tion: 


Will  the  LQG/LTR  technique  of  controller  design  produce  an  effective 
controller  for  a  small  satellite  using  a  non-traditional  elastic  memory 
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composite  flexible  appendage  and  how  robust  is  the  system  in  the  face  of 
material  uncertainties,  high  frequency  modeling  errors,  and  appendage 
length  modifications? 

1.6  Research  Importance 

This  research  is  both  innovative  and  significant.  A  brief  summary  of  the  literature 
search  in  the  area  of  flexible  spacecraft  modeling  is  provided  here  while  a  detailed 
literature  study  is  provided  in  earlier  sections  within  this  chapter.  This  brief  summary 
is  given  to  set  the  stage  for  the  reader  on  what  researchers  have  done  in  the  area 
of  flexible  spacecraft  modeling.  References  are  provided  in  author,  year  format  while 
specific  reference  numbering  is  provided  in  the  detailed  literature  search  discussion. 

1.6.1  Innovative 

Modeling  flexible  spacecraft  is  not  a  new  area  of  study.  60  years  ago,  structural 
analysis  studies  were  first  applied  to  spacecraft  at  the  beginning  of  the  ”  space  race” 
between  the  United  States  and  the  former  Soviet  Union.  Likins  developed  hybrid  coor¬ 
dinate  equations  using  cantilever  boundary  conditions  in  1970,  which  was  incorporated 
into  NASA’s  Space  Vehicle  Design  Criteria  documents  in  1971.  Cantilever  boundary 
conditions  were  continued  to  be  utilized  in  the  papers  surveyed  by  Croopnick  et  al.  in 
1979,  by  Meirovitch  in  1979,  Balas  in  1982,  and  Nurre  et  al.  in  1984.  The  focus  of 
the  papers  included  in  all  of  the  surveys  focused  on  the  control  of  large  flexible  space 
structures  and  used  cantilever  boundary  conditions  in  the  formulation  of  the  system’s 
equations  of  motion.  Some  researchers  mentioned  that  a  similar  process  can  be  done  for 
free-free  boundary  conditions,  but  none  of  them  demonstrated  or  implemented  equa¬ 
tions  of  motion  using  free-free  mode  shapes;  Canavin  1977,  Sundararajan  1987,  Junkins 
1993,  and  Izzo  2004. 

Over  the  years,  spacecraft  have  become  smaller  in  design.  Even  though  several 
small  flexible  spacecraft  would  best  be  characterized  as  free-free  systems,  cantilever 


29 


assumptions  are  still  being  used  (Gorinevsky  1998,  Maghami  1998,  Wie  1998,  Lomas 
2001,  Bodineau  2004,  Calise  2004,  Gili  2004,  Izzu  2004,  Tafazoli  2004,  Takahito  2004, 
Wilson  2004,  and  Ledesma  2005).  This  researcher  is  not  stating  these  well  informed 
people  are  making  large  errors  in  modeling.  In  some  of  the  cases  presented,  assuming 
cantilever  mode  shapes  is  a  valid  assumption,  but  this  assumption  cannot  be  used  on  the 
formulation  of  the  equations  of  motion  for  all  small  flexible  spacecraft  (see  Section  2.2.2). 

This  research  applies  free-free  boundary  conditions  to  a  small  flexible  space¬ 
craft  using  a  non-traditional  gravity  gradient  boom  constructed  from  Elastic  Memory 
Composite  material,  demonstrates  the  development  of  the  equations  of  motion  using 
these  conditions,  and  implements  the  equations  in  an  attitude  control  scenario  using 
LQG/LTR  techniques.  A  key  finding  is  the  historical  assumption  of  considering  the 
spacecraft  as  a  rigid  body  if  the  first  resonant  modes  are  greater  than  an  order  of  mag¬ 
nitude  when  compared  to  the  controller  bandwidth  may  no  longer  be  a  valid  assumption 
when  a  small  flexible  spacecraft  is  taken  into  consideration. 

Two  methods  are  described  and  compared  for  reducing  the  order  of  the  con¬ 
trollers.  Reducing  the  order  of  the  controllers  leads  to  fewer  coefficients  needed  to  be 
coded  onto  the  satellite’s  on-board  processor.  Comparisons  of  the  number  of  opera¬ 
tions  per  computing  cycle  are  discussed  to  aid  a  satellite  designer  during  the  on-board 
processor  sizing  task.  In  addition,  several  tools  are  provided  for  small  satellite  designers 
to  determine  the  affects  uncertainties  in  the  material  as  well  as  varying  the  length  of 
a  gravity  gradient  boom  will  have  on  the  satellite’s  performance  parameters  (such  as 
settling  time  and  stability  robustness  values). 

1.6.2  Significance 

The  intention  of  this  research  is  to  provide  an  understanding  of  how  attitude  con¬ 
trol  of  small  satellites  is  impacted  when  attached  flexible  appendages  are  constructed 


from  materials  which  are  more  elastic  than  traditional  structures.  A  better  under- 
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standing  of  how  variations  in  elasticity  and  length  of  these  appendages  impact  attitude 
control  authority  will  aid  designers  in  reducing  some  of  the  risks  involved  in  using  soft 
materials. 

Although  the  2002  CTD  study  makes  EMC  materials  look  promising,  there  are 
concerns.  To  date,  only  deployment  hinges  have  been  tested  by  CTD.  In  2004,  CTD 
developed  a  few  test  longerons  constructed  from  EMC  materials.  However,  only  ther¬ 
mal  cycling  tests  of  small  lengths  (approximately  12  inches)  in  a  vacuum  chamber  were 
conducted  to  determine  battery  sizing  requirements  to  heat  the  material  prior  to  deploy¬ 
ment.  Initial  dynamic  analysis  of  an  EMC  appendage  offers  preliminary  data;  however, 
a  system  test  of  an  engineering  model  of  the  deployable  structure  has  not  been  con¬ 
ducted.  This  analysis  indicates  a  deployable  structure  made  from  EMC  materials  will 
have  a  bending  frequency  of  1.5  Hz  and  a  torsional  frequency  of  1.7  Hz. 

Constructing  a  deployable  appendage  from  EMC  material  is  a  novel  idea.  The 
deployment  mechanism  is  contained  within  the  stored  stress  energy  of  the  material  itself 
and  does  not  require  mechanically  complicated  motors  or  added  inhibit  precautions 
needed  on  currently  used  traditional  style  booms.  Small  satellites  will  take  advantage 
of  the  lighter  mass  systems  and  will  begin  to  utilize  missions  requiring  longer  boom 
lengths. 

The  proposed  research  will  support  both  commercial  and  governmental  areas  of 
research  and  development.  Demonstrating  the  feasibility  of  space  structures  constructed 
from  EMC  materials  is  the  next  step  in  CTD’s  current  design  and  implementation  ef¬ 
forts  in  the  area  of  EMC  space  structures.  The  Space  Vehicles  Directorate  of  the  Air 
Force  Research  Labs  is  interested  in  CTD’s  technological  advancements.  In  particular, 
the  Power  Sail  program  is  considering  the  use  of  EMC  deployment  systems  in  their  con¬ 
ceptual  design.  This  study  will  identify  the  key  areas  of  effectiveness  in  controlling  such 
flexible  space  structures  which  will  prove  vital  in  allowing  AFRL  to  pursue  advanced 
designs  in  their  Power  Sail  structure. 
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The  United  States  Air  Force  Academy  will  use  the  findings  of  this  research  to 
complete  risk  assessments  in  the  area  of  attitude  determination  and  control  of  an  EMC 
gravity  gradient  boom  providing  passive  control  to  a  small  satellite.  A  successful  in¬ 
tegration  of  the  boom  into  the  next  FalconSat  spacecraft  may  provide  an  opportunity 
to  generate  flight  heritage  with  the  system  before  the  end  of  2010.  This  will  offer  a 
new  way  of  deploying  flexible  space  structures  without  the  disadvantages  of  traditional 
systems. 


Chapter  2 


Equations  of  Motion 


2.1  Flexible  Spacecraft  Modeling 

The  attitude  motion  of  a  flexible  spacecraft  is  properly  described  by  coupled  sets 
of  partial  and  ordinary  differential  equations.  The  rotational  motion  of  the  undeformed 
system,  called  the  rigid  body  motion,  is  described  by  ordinary  differential  equations 
while  the  flexures  are  described  by  partial  differential  equations.  The  rigid  body  dy¬ 
namics  are  derived  from  Euler’s  rotational  equations  of  motion  and  include  gravity 
gradient  torques.  Numerical  finite  element  models  are  used  to  determine  the  mode 
shapes  and  natural  frequencies  as  well  as  the  mass,  damping,  and  stiffness  matrices  of 
the  flexible  system [30].  The  assumed  modes  method  is  used  to  couple  the  rigid  body 
and  flexible  dynamics  by  using  the  spatial  solutions  of  the  partial  differential  equations 
as  assumed  mode  shapes  and  letting  the  modal  coefficients  serve  as  the  generalized 
coordinates  describing  the  flexures. 

This  research  considers  the  impact  a  flexible  appendage  has  on  the  attitude  control 
system  of  a  small  satellite.  Topics  such  as  meeting  pointing  requirements  and  attitude 
maneuvers  are  concerned  with  the  rotational  motion  of  the  spacecraft  while  station 
keeping  and  changes  in  orbital  parameters  deal  with  the  translation  of  the  spacecraft. 
Eq.  A. 62  from  Appendix  A  illustrates  how  the  translational  effects  can  be  removed  from 
the  vibrational  equations  to  leave  a  coupled  set  of  rotation-vibration  equations. 
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2.1.1  Rigid  Body  Dynamics 

To  begin  the  formulation  of  the  equations  of  motion  for  a  small  satellite  with  a 
flexible  appendage,  consider  Euler’s  rotational  equations  of  motion  shown  as  Eq.  C.34 
in  Appendix  C. 

J\iOi  —  ( J-2  —  J3)^2^3  =  Td 
J2CO2  —  (:h  —  -h  Vi^’3  =  TC2 

J3U3  ~  (’h  —  J2V1W2  =  Tc 3  (2.1) 

where  J2,  and  J3  are  the  principal  moments  of  inertia  of  the  undeformed  system, 
and  lo3  are  the  angular  rates  of  motion  about  the  principal  axes  (yaw,  pitch,  and 
roll  respectively),  and  Tci,  Tc 2,  and  Tc 3  are  the  attitude  control  torques  of  the  spacecraft 
about  these  axes. 

Expanding  out  the  coupled  terms,  LOiLOj,  and  writing  Eq.  2.1  in  matrix  notation 
produces 
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V ic3  J 


Since  the  spacecraft  being  modeled  is  on  orbit  and  the  flexible  appendage  is 
providing  passive  gravity  gradient  stabilization,  gravitational  forces  need  to  be  included. 
Using  the  local  vertical  and  local  horizontal  (LVLH)  reference  frame  to  describe  the 
orientation  of  the  spacecraft  places  the  1st  axis  along  the  orbit  direction  (velocity  or 
ram  direction),  the  2nd  axis  perpendicular  to  the  orbital  plane  (orbit  normal  direction), 
and  the  3rd  axis  pointing  towards  the  Earth  (nadir  pointing). 

Using  the  gravity  gradient  derivations  included  in  Appendix  B  and  the  rotational 
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kinematics  shown  in  Appendix  C,  the  right  hand  side  of  Eq.  2.2  become 
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where  n  = 


+  3n2  |  (J,  -  JZ)C^CW  I  (2-3) 

^  (— +  -hC\:iC2s  J 

is  the  orbital  rate,  /i  =  3.9860xl014m3/sec2  is  the  gravitational  para¬ 
meter  for  Earth,  and  R0  is  the  orbital  radius  of  the  spacecraft  measured  from  the  center 
of  the  Earth  to  the  center  of  mass  of  the  spacecraft.  Ci3,C2z  and  C33  are  elements  of 
the  transformation  matrix  used  to  go  from  the  LVLH  frame  to  the  body  fixed  frame. 

Writing  the  equations  of  motion  in  differential  form  and  moving  the  gravity  gra¬ 
dient  torques  to  the  left  hand  side  yields: 


Jld>l  —  (J2  —  J3)u)2(jJ3  +  3n2(J2  —  J3)C2sC33  =  Tc  1 

,J2uj2  —  (J3  —  J\)ujiuj3  +  3n2(,/3  —  .J\)C\  3  C33  =  Tc  2 

^3^3  —  (Ji  —  J2)wilu2  +  3?r2(Ji  —  J2)Ci3C23  =  Tc  3 


(2.4) 


As  shown  in  Appendix  C.1.2,  a  singularity  occurs  in  the  kinematic  equations.  For 
a  gravity  gradient  stabilized  satellite,  it  is  assumed  the  boom  will  be  either  zenith  or 
nadir  pointing.  The  only  way  the  pitch  or  roll  angles  can  equal  90°  is  if  the  boom  is 
perpendicular  to  the  gravitational  force  vector. 

Consider  the  rotational  sequence  of  C\(B\)  <—  C2(B2)  <—  03(63)  to  the  body  frame 
from  the  LVLH  frame: 
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The  angular  velocity  for  this  rotational  sequence  is 
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To  solve  for  the  kinematic  differential  equations,  the  3x3  non-orthogonal  matrix 
above  is  inverted  to  finally  yield: 
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where  the  singularity  now  occurs  with  a  pitch  value  of  90°. 

Going  back  to  the  equations  of  motion  in  Eq.  2.4  and  inserting  the  values  C'1,3  = 
—sO 2,  C23  =  s9ic02,  and  C33  =  c0ic02  from  Eq.  2.5,  the  equations  of  motion  become 

J\Co\  —  (J2  —  Jz)u)2U3  +  3n2(J2  J‘s)s0\C02 C0\C02  =  Tc  1 

J2W2  —  (J3  —  Ji)wicn3  —  3?r2(J3  —  Ji)s02C0ic02  =  TC2 

<^3^3  —  (J\  —  J2)^lL02  —  3n2(Ji  —  J2)s02S0iC02  =  TC3  (2-13) 

Recall  that  the  intention  of  this  research  is  to  determine  the  effect  the  flexible 

appendage  has  on  the  attitude  control  authority  of  a  small  satellite.  Large  slewing 

maneuvers  of  the  satellite  are  not  considered  while  certain  pointing  requirements  are 
maintained.  Therefore,  it  is  safe  to  assume  the  satellite  is  moving  through  small  angular 
displacements. 

s9  «  9 

c9  «  1 

0X02  ~  0 


Applying  small  angle  approximation  to  Eq.  2.13  yield: 

Ji^i  ~  (J2  ~  J'3)^2^3  +  3n2(J2  —  Jz)9\  =  Tc  1 
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Applying  small  angle  approximation  to  the  kinematic  equations,  Eq.  2.11,  yield: 
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Eq.  2.15  in  differential  form  is 
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u)2  =  0-2  +  6*16*3  -  n 

w3  =  —O1O2  +  #3  +  nO  1 

(2.16) 

but  OiOj  «  0,  so 


OJl 

=  6\  —  n6*3 

io2 

=  02  -  n 

CJ3 

=  93  +  n6*i 

(2.17) 

If  orbital  motion  is  assumed  to  be  constant  (i.e.  no  eccentricity  in  orbit),  the  time 
differential  of  Eq.  2.17  is 

uj\  =  9 1  —  n#3 

U>  2  =  92 

(j3  =  03  +  n6i  (2.18) 

Inserting  Eq.  2.17  and  Eq.  2.18  into  Eq.  2.14,  applying  small  angle  approximation, 

and  collecting  like  terms  results  in 

+  n(-Ji  +  J2-  h%  +  An\j2-  J3)0i  =  Tcl 
J202  -  37i2(J3  -  Ji)02  =  Tc2 

^3^3  +  —  J2  +  ^3)^1  —  n2{Ji  ~  ^2)^3  =  Ted,  (2-19) 

Eq.  2.19  is  the  rigid  body  dynamics  of  a  satellite  in  a  circular  orbit  with  gravity 
torques  included.  The  dynamics  of  the  flexible  appendage  still  need  to  be  included. 

2.1.2  Flexible  Dynamics 

The  method  used  to  include  the  flexible  dynamics  of  the  appendage  with  the  rigid 
body  of  motion  of  the  satellite  transforms  the  equations  of  motion  in  physical  coordinates 
to  decoupled  vibrational  equations  (also  known  as  modal  equations  but  renamed  here  to 
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prevent  confusion  with  the  modal  form  of  the  state  space  equations) .  This  is  done  using 
a  linear  coordinate  transformation  known  as  the  modal  transformation.  In  Appendix  A, 
the  displacement  of  the  flexible  system,  u(r,t),  is  expressed  as 

OO 

u(r,t)  =  <t>m(r)r}m(t)  (2.20) 

m= 0 

where  rjm(t )  is  the  mth  vibrational  coordinate  and  (j>m(r)  is  the  mth  normal  mode  shape 
of  the  mode  shape  matrix,  [</>],  whose  columns  are  the  eigenvectors  of  the  system. 

The  differential  equations  of  motion  for  an  undamped  free  vibration  system  are 
of  the  form 

[m]il  +  [k]u  =  0  (2-21) 

where  the  matrices  [m]  and  [A;]  are  arbitrary  mass  and  stiffness  matrices  of  the  system 
with  symmetric  and  constant  elements.  Since  [<; t>\  is  a  spatial  variable  and  rj  is  a  temporal 
variable,  we  can  insert  Eq.  2.20  into  Eq.  2.21  and  get 

[m]  [cj)]ii  +  [k]  [4>\r)  =  0  (2.22) 

Premultiplying  by  [<j)]T  yields 

[M}t)  +  [Kjr]  =  0  (2.23) 

where 

[A-/]  =  =  [M]T 

m  =  wT\m  =  wT 

If  [(/>]  is  orthonormal,  then  the  generalized  mass  matrix,  [M\,  is  the  identity  ma¬ 
trix  ([M]  =  diag(l.  •••,!))  and  the  generalized  stiffness  matrix,  [K],  is  a  diagonal  matrix 
whose  elements  are  equal  to  the  natural  frequencies  squared  ([A"]  =  diag{uQ,uj\,  ■  ■  ■ ,  u>^)). 
Since  both  generalized  matrices  are  diagonal,  it  makes  the  vibrational  equations  both 
inertially  uncoupled  and  elastically  uncoupled.  This  results  in  independent  equations 
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of  flexible  motion.  If  there  are  m  mode  shapes  included  in  the  analysis,  these  equations 
look  like  the  following: 

r'n  +  ufm  =  0 

m  +  u\r)2  =  0 

Vm  T  W mVm  =  0  (2.24) 

where  r]m  is  the  mth  vibrational  coordinate  associated  with  the  mth  mode  shape  and 
ojm  is  the  mth  natural  frequency.  Therefore,  r ]  has  dimensions  of  (lxm).  The  mode 
shape  matrix,  [(j)\  will  have  dimensions  of  (6nxm)  where  n  is  the  number  of  nodes 
used  in  the  finite  element  analysis  and  each  node  has  six  degrees  of  freedom  (three 
translation  and  three  rotation  for  each  node).  The  number  of  nodal  points  used  in 
the  FEM  analysis  is  up  to  the  discretion  of  the  researcher.  However,  a  more  accurate 
estimate  of  the  deformation  of  the  flexible  system  is  found  by  using  more  nodal  points 
to  reduce  approximation  errors  in  the  displacement  of  each  element. 

The  vibrational  equations  of  motion  are  propagated  along  with  the  rigid  body 
equations.  If  the  movement  of  certain  points  along  the  flexible  appendage  are  of  in¬ 
terest,  then  inserting  the  updated  vibrational  coordinates  into  Eq.  2.20  will  generate 
the  displacement  function  u(r ,  t )  in  physical  coordinates.  To  illustrate  this,  consider  a 
simple  example  of  a  flexible  appendage  in  which  three  nodes  are  selected  for  the  FEM 
analysis. 

Each  node  is  free  to  translate  in  three  directions  (Tl,  T2,  T3)  but  for  now  we 
will  disregard  rotation  in  three  directions  (Rl,  R2,  R3).  Also,  let’s  assume  three  mode 
shapes  are  used  in  the  analysis.  Then  the  physical  coordinates  of  the  appendage  are 
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found  as  follows: 
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where,  for  example,  Node3TlMS2  is  the  eigenvector  value  describing  the  displacement 
of  the  third  node  in  the  first  direction  in  response  to  the  second  mode  shape.  If  the 
total  displacement  of  the  third  node  in  the  first  direction  is  required,  it  would  be 


U  =  <h,lVl  +  07,2m  +  07,3113  (2-26) 

This  can  be  done  for  each  of  the  three  nodes  in  the  three  translational  directions. 
More  complicated  models  which  contain  a  greater  number  of  nodes  will  produce  better 
estimates  of  the  displacement  shape  of  the  flexible  system,  but  will  require  more  nodes 
as  well  as  rotation  in  three  directions  along  with  translational  motion. 


2.1.3  The  Coupling  Equation 

The  rigid  body  equations  need  to  be  coupled  with  the  vibrational  equations. 
The  hybrid  coordinate  approach[129]  presents  the  rigid  body  equations  and  appendage 
deformation  in  the  form 

Id-dTi]  =  Tc 
ij  +  2( u:nf]  +  ufy ?  -  86  = 


0 


(2.27) 
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where  5  is  called  the  coupling  matrix[127].  Solving  the  Lagrange  equations  in  Appen¬ 
dix  A  yields  the  following  definition  for  the  coupling  matrix: 


5  =  Xf  R  —  Xj 


(2.28) 


where 


A,  = 

A3  = 


/  app 


[4>\  dm  (3xm) 
r  [(f)]  dm  (3  x  m) 


(2.29) 


>  app 


Closed  form  solutions  can  be  found  for  the  mass  and  inertia  integrals  in  Eq.  2.29 
and  the  coupling  matrix  is  rewritten  as 


5  =  —rfPMc^oE  ~  E eoR  ~  rY>Eo) 


(2.30) 


where  (/>  is  the  truncated  mode  shape  matrix,  Mc  is  the  generalized  (6Nx6N)  inertia 
matrix  of  N  cantilevered  appendages,  T,oe  and  T^eo  are  summation  matrices  consisting 
of  ones  and  zeros,  and  R  and  f  are  skew  symmetric  matrices  of  position  vectors.  Each 
variable  is  explained  in  further  detail  below  to  aid  in  the  application  of  Eq.  2.30. 

A  concern  arises  when  looking  at  the  dimensions  of  [ (j>\  used  in  Eq.  2.29.  As  shown, 
the  coupling  matrix  has  dimensions  (mx3).  However,  the  mode  shape  matrix  takes  the 
form  (6nxm).  The  coupling  matrix  couples  the  rotational  motion  of  the  satellite  with 
the  flexible  motion  of  the  appendage  and  there  is  a  constraint  relationship  between  the 
two  as  illustrated  in  Figure  2.1. 

Ug/c  is  a  point  located  on  the  spacecraft  while  Uapp i  is  the  location  of  the  first 
node  used  in  the  FEM  analysis.  Notice  that  both  of  these  points  are  located  at  the 
connection  point  between  the  spacecraft  and  the  flexible  appendage.  As  a  result,  the 
displacement  and  rotation  of  Ug/c  is  the  same  as  those  of  Uapp i-  It  is  at  this  location 
where  rotational  (and  translational)  energy  is  passed  between  the  spacecraft  and  the 
flexible  appendage.  Therefore,  while  the  mode  shape  matrix  may  take  the  form  of 


Spacecraft 


Figure  2.1:  Rotational- Vibrational  Constraint  Illustration 
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the  coupling  matrix  only  relies  on  the  translations  and  rotations  of  the  first  node. 
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which  now  puts  the  mode  shape  matrix  in  a  truncated  form  with  dimensions  (6xm). 
Note,  Eq.  2.32  contains  information  for  only  one  appendage.  If  N  additional  flexible 
appendages  are  considered,  append  the  rows  of  the  truncated  mode  shape  matrix  with 
the  nodal  information  for  the  connection  points  to  each  flexible  structure.  Thus,  the 
dimensions  become  (6Nxm). 

The  generalized  inertia  matrix  of  cantilevered  appendages,  Mc,  is  comprised  of 
the  masses  and  inertias  of  each  connected  appendage  and  takes  the  form 


/ 


Mr  = 


ml  0000  0  0 

0  /1  0  0  0  0  0 

0  0  m2  0  0  0  0 

0  0  0  /2  0  0  0 

0  0  0  0  0  0 


\ 


0 


0  0  0  0  rriN  0 


V 


0  0  0  0  0 


0  IN 


(2.33) 


where  ml  is  the  (3x3)  mass  matrix  of  the  first  appendage  such  that 

/  -  \ 


ml  = 


rrhipp  1  0  0 

0  m.appi  0 

^  0  0  mappi  J 


(2.34) 


44 


and  II  is  the  (3x3)  inertia  matrix  of  the  first  appendage  such  that 


/ 


II  = 


I appx  0 


V 


appy 


0 


0 

0 


\ 


(2.35) 


*■  appz  J 

R  and  r  are  skew  symmetric  matrices  of  the  position  vectors  which  have  the 
general  form 


/ 


Rl  = 


\ 


V 


(2.36) 
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o  -ri3  ri2 

Rl3  0  -Rli 
-Rl2  Rl  i  0 

where  Rl  is  the  position  vector  from  the  system  center  of  mass  to  to  the  connection 
point  of  the  first  appendage  and  rl  is  the  position  vector  from  the  connection  point  of 
the  first  appendage  to  the  center  of  mass  of  the  first  appendage. 

Multiple  appendages  are  accounted  for  by  utilizing  the  summation  matrices  Hoe 
and  Heo ■  These  matrices  are  Boolean  operator  matrices  of  the  form 

T 


Hoe  = 
Heo  = 


O  E  O  E 
E  O  E  O 


O  E 


E  O 


(2.37) 


where  O  is  a  (3x3)  zero  matrix  and  E  is  a  (3x3)  identity  matrix. 

The  coupling  matrix  is  general  enough  that  the  same  concept  not  only  applies  to 
a  satellite  with  a  flexible  gravity  gradient  boom,  but  also  works  for  solar  arrays,  multiple 
flexible  appendages  of  various  shapes  and  elasticity,  etc.  What  changes  from  each  of 
these  instances  is  the  mode  shape  matrix  generated  during  FEM  analysis.  The  coupling 
matrix  contains  all  of  the  dynamics  used  in  coupling  the  rigid  body  dynamics  with  each 
of  the  flexible  appendages. 


Table  2.1:  Material  Properties  of  Beryllium  Copper  and  EMC 
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Beryllium  Copper 

EMC 

Young’s  Modulus,  E  [GPa] 

138 

40.06 

Poisson  Ratio,  v 

0.30 

0.31 

Shear  Modulus,  G  [GPa] 

53.1 

15.29 

Density,  p  [||] 

8830 

1384 

2.2  Finite  Element  Model 

Since  the  coupling  matrix  is  dependent  upon  the  properties  of  the  satellite,  the 
type  and  number  of  flexible  appendages,  and  the  materials  used  in  the  construction  of 
the  appendage,  the  illustration  shown  in  Figure  2.1  will  describe  the  set  up  of  the  system 
to  be  modeled.  The  spacecraft  will  be  a  typical  small  satellite  based  on  the  FalconSat 
system  described  in  Appendix  D.  The  flexible  appendage  is  a  non-traditional  gravity 
gradient  boom  constructed  from  Elastic  Memory  Composite  (EMC)  materials  with  a 
cubic  tip  mass. 

2.2.1  Model  Parameters 

As  mentioned  in  Section  1.2,  beryllium  copper  is  a  common  alloy  used  in  tra¬ 
ditional  beam  elements.  Table  2.1  lists  the  material  properties  for  both  beryllium 
copper [1]  and  EMC  materials [206].  The  composite  material  used  in  the  finite  element 
model  (FEM)  not  only  is  more  flexible  but  also  has  15.7%  the  density  of  materials  which 
traditional  booms  are  constructed  from. 

The  EMC  material  properties  were  entered  in  MSC.Patran  2004  r2  to  create  a 
FEM  of  a  cantilevered  beam  with  lumped  masses  at  either  end  to  represent  the  satellite 
and  tip  mass.  The  satellite  lumped  mass  boundary  conditions  are  fixed  while  the  tip 
mass  is  allowed  to  move.  The  lumped  masses  were  treated  as  cubic  shapes  with  uniform 
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Table  2.2:  Lumped  Mass  Properties  Used  in  FEM  Analysis 


Satellite 

Tip  Mass 

Mass 

40% 

7kg 

Side  Length 

0.5m. 

0.25  m 

Inertia 

1.667%m2 

0.0729 kgm'2 

mass  distribution.  For  a  cubic  shape  with  uniform  mass  distribution,  the  inertia  about 
the  object’s  center  of  mass  is 

hi  =  I22  =  h:i  =  —  (2.38) 

6 

where  m.  is  the  mass  of  the  object,  l  is  the  length  of  a  side  and  the  products  of  inertia 
are  equal  to  zero.  Table  2.2  lists  the  properties  entered  into  the  FEM  for  the  satellite 
and  tip  mass. 

The  beam  representing  the  EMC  flexible  appendage  is  modeled  as  a  hollow,  cylin¬ 
drical  tube  with  properties  shown  in  Table  2.3  to  resemble  the  initial  design  character¬ 
istics  of  the  EMC  appendage [206].  One  hundred  node  points  were  used  in  the  analysis 
evenly  spaced  4cm  apart  to  minimize  the  approximation  errors  in  the  mode  shapes. 

A  normal  modes  solution  type  was  selected  to  complete  a  full  run  of  the  entire 
model.  This  solution  will  generate  the  eigenvectors  for  all  available  modes.  The  nominal 
(or  truth)  model  used  in  the  analysis  includes  the  first  eighteen  resonant  modes,  regard¬ 
less  of  their  impact  (as  with  the  compression  mode)  or  modes  which  are  disregarded 
during  order  reduction  methods  of  the  controller. 

The  first  eighteen  mode  shapes  generated  for  the  cantilevered  appendage  are  listed 
in  Table  2.4  while  plots  of  the  first  three  bending  modes  and  the  torsional  mode  are 
provided  in  Figures  2. 2-2. 5.  One  thing  to  consider  is  the  compression  mode  places  very 
minimal  torque  on  the  satellite  if  the  center  of  mass  of  the  appendage  is  placed  on 
one  of  the  principal  axes  of  the  system.  This  is  validated  by  looking  at  the  angular 
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Table  2.3:  Beam  Properties  Used  in  FEM  Analysis 


Property 

Value 

Total  Length 

4m 

Outer  Diameter 

2.54  cm 

Thickness 

0.46482 mm 

Cross  Sectional  Area 

73.5  mm2 

hi  =  I'22 

2.3  x  10 ~6kgm2 

displacement  eigenvector  of  the  compression  mode  for  the  first  node  of  the  appendage 
connected  to  the  satellite. 

/  - - \  (  2.446  x  10"15  ^ 

2.522  x  10"13 
-2.765  x  10-16 


V 


/ 


(2.39) 


N1R1MS8 
N1R2MS8 
N1R3MS8 

The  compression  mode  is  left  in  the  nominal  model,  but  will  most  likely  not 


V 


/ 


appear  in  the  controller  design  model  because  of  order  reduction  techniques. 
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Table  2.4:  Cantilever  Resonant  Frequencies  From  FEM  Analysis 


Mode  Shape 

Frequency  (Hz) 

1st  Bending  in  Y  direction 

0.3947 

1st  Bending  in  X  direction 

0.3947 

1st  Torsional  Mode 

7.8585 

2nd  Bending  in  Y  direction 

12.141 

2nd  Bending  in  X  direction 

12.141 

3rd  Bending  in  Y  direction 

29.322 

3rd  Bending  in  X  direction 

29.322 

1st  Compression  Mode 

51.116 

4th  Bending  in  Y  direction 

61.567 

4th  Bending  in  X  direction 

61.567 

5th  Bending  in  Y  direction 

115.62 

5th  Bending  in  X  direction 

115.62 

6th  Bending  in  Y  direction 

188.56 

6th  Bending  in  X  direction 

188.56 

7th  Bending  in  Y  direction 

279.26 

7th  Bending  in  X  direction 

279.26 

8th  Bending  in  Y  direction 

386.96 

8tfl  Bending  in  X  direction 

386.96 
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Figure  2.2:  First  Bending  Mode  For  Cantilever  Conditions 
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Figure  2.3:  Second  Bending  Mode  For  Cantilever  Conditions 
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Figure  2.4:  Third  Bending  Mode  For  Cantilever  Conditions 
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Fringe:Default  A1  :Mode  3  :  Freq.  =  7.8585:  Eigenvectors,  Translational-(NON-LAVERED)  (MAG) 
Deform  Default  A1  :Mode  3  :  Freq.  =  7.8585:  Eigenvectors,  Translational 
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Figure  2.5:  First  Torsional  Mode  For  Cantilever  Conditions 
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Figure  2.6:  Cantilever  Beam  with  Applied  Point  Load 


2.2.2  Free-Free  Boundary  Conditions 

The  traditional  hybrid  coordinate  equation  model  first  presented  by  Likins  in  a 
1970  Jet  Propulsion  Laboratory  Technical  Report [127]  has  strongly  influenced  the  dy¬ 
namical  modeling  of  flexible  spacecraft.  The  simplified  rotational/vibrational  coupled 
equations,  Eq.  2.27,  is  found  in  one  form  or  another  in  the  works  of  several  key  re¬ 
searchers  in  the  field  of  flexible  dynamics.  While  the  references  are  too  numerous  to  list 
here,  the  survey  papers  presented  in  [177]  [50]  [149]  [15]  [161]  provide  a  broad  enough 
summary  of  the  collection  of  work. 

The  literature  considers  cantilever  modes  when  analyzing  the  flexible  appendage 
(see  Section  1.6.1).  A  common  definition  in  the  structures  community  for  a  cantilever 
beam  is  a  slender  member  which  is  fixed  at  one  end  and  free  at  the  other  that  supports 
loadings  that  are  applied  perpendicular  to  their  longitudinal  axis [86].  The  boundary 
conditions  applied  to  the  fixed  end  of  the  beam  restrain  both  translational  and  rotational 
motion.  During  the  generation  of  the  equations  of  motion  for  the  beam,  either  a  point 
load  is  applied  at  some  location  along  the  appendage,  Figure  2.6,  or  as  a  distributed 
load  along  the  appendage  length,  Figure  2.7. 

When  a  flexible  appendage  is  attached  to  a  symmetrical  spacecraft,  the  root 
node,  or  connection  point  between  the  appendage  and  spacecraft  bus,  does  both  ro¬ 
tate  and  translate.  If  one  looks  at  the  system  as  a  whole,  cantilever  modes  doesn’t 
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Figure  2.7:  Cantilever  Beam  with  Applied  Distributed  Load 


make  sense  due  to  the  motion  of  the  root  node.  The  literature  assumes  the  flexible 
appendage  is  ’’cantilever  like”  when  deriving  the  mode  shape  matrix  and  the  accepted 
practice  in  flexible  spacecraft  modeling  is  to  calculate  the  cantilever  mode  shapes  of  the 
attached  appendage  and  couple  the  resulting  mode  shape  matrix  with  the  spacecraft 
non-vibrational  equations  of  motion.  This  hybrid  approach  is  valid  if  one  looks  at  a 
body  fixed  reference  frame  located  at  the  root  node,  since  in  this  frame  the  connection 
point  between  the  appendage  and  spacecraft  bus  does  not  appear  to  either  rotate  or 
translate. 

Modeling  the  flexible  appendage  as  a  cantilevered  beam  is  the  correct  approach 
when  dealing  with  attitude  control  of  large  space  structures.  The  boundary  conditions  at 
one  end  have  the  appendage  fixed  in  translation  and  rotation.  This  is  a  valid  assumption 
when  the  total  system  center  of  mass  is  located  close  to  the  center  of  mass  of  the 
controlling  body.  Two  examples  of  when  this  occurs  is  large  space  structures  when 
mcB  mapp  and  when  flexible  appendages  are  symmetrically  orientated  about  the 
controlling  body  (as  shown  in  Figure  2.8  and  Figure  2.9).  The  torque  generated  by  the 
displaced  appendage  is  applied  at  the  connection  point  between  the  controlling  body 
and  the  appendage.  The  torque  experienced  by  the  controlling  body  is  in  the  same 
direction  as  it  is  generated  from  the  appendage. 

For  both  the  large  flexible  spacecraft  and  symmetrical  appendage  configuration 
examples,  the  position  of  the  system  center  of  mass  experiences  small  displacements [31]. 
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Figure  2.8:  Three-Axis-Stabilized  Geosynchronous  Communications  Satellite[237] 


Figure  2.9:  Example  of  a  Symmetrical  Appendage  Configuration 
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As  the  system  is  deformed  in  asymmetrical  modes,  the  velocity  of  the  system  center  of 
mass,  and  therefore  the  time  rate  of  change  of  the  moments  of  inertia,  are  approximately 
zero  [100].  This  eliminates  several  terms  from  the  Lagranian  equations  and  simplifies  the 
resultant  equations  of  motion.  Likins’  work  retained  asymmetric  modes  corresponding 
to  known  cantilever  modes  while  showing  that  the  symmetric  modes  contributed  nothing 
and  could  be  ignored[129]. 

For  a  gravity  gradient  stabilized  small  satellite,  the  assumption  of  a  cantilever 
appendage  is  no  longer  as  valid  as  it  is  for  a  large  space  structure.  The  flexible  appendage 
is  attached  to  both  the  spacecraft  bus  and  the  tip  mass.  To  determine  which  end  is 
free,  one  can  argue  that  the  end  which  does  not  generate  loads  perpendicular  to  the 
longitudinal  axis  would  be  free  to  move.  Thus,  the  tip  mass  is  merely  considered  as  part 
of  the  appendage  which  is  no  longer  uniform  in  mass  distribution.  Tethered  spacecraft, 
however,  confuse  this  definition  if  both  masses  tethered  are  identical  each  with  their 
own  attitude  control  systems.  In  these  cases,  it  is  no  longer  clear  which  end  is  fixed  and 
which  end  is  free.  The  conventional  definition  of  cantilevered  beams  breaks  down. 

This  is  especially  true  as  the  mass  of  the  satellite  gets  closer  to  the  mass  of  the 
attached  appendage  and  tip  mass.  As  these  two  ends  come  closer  to  each  other  in  mass, 
the  system  center  of  mass  moves  away  from  a  location  within  the  constrained  boundary 
condition  and  is  found  at  some  point  along  the  appendage  itself.  A  better  approach  is  to 
consider  the  appendage  as  a  free-free  beam  with  attached  lumped  masses  at  either  end 
for  the  satellite  and  tip  mass.  Now,  each  end  of  the  appendage  is  allowed  to  translate 
and  rotate  when  determining  the  mode  shapes  used  to  couple  the  rigid  body  and  flexible 
dynamics.  This  is  a  more  realistic  assumption  when  dealing  with  small  satellites  because 
the  attitude  of  the  satellite  is  directly  impacted  by  the  motion  of  the  appendage.  For  a 
cantilever  assumption,  the  satellite  doesn’t  move  as  a  result  of  appendage  motion. 

For  the  free-free  assumption,  the  center  of  mass  of  the  total  system  is  not  nec¬ 
essarily  contained  within  the  controlling  body,  nor  along  one  of  the  principle  axes.  In 
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Table  2.5:  Free-Free  System  Eigenvalues  Using  Unmodified  HCEs 


Eigenvalue 

Coupled  Pair 

0 

0 

-0.19803  +  19.8023j 

-0.19803  -  19.8023j 

0  +  0.00217j 

0  -  0.00217j 

0  +  0.00188j 

0  -  0.00188j 

67.0388 

-67.0388 

9.08888 

-9.08888 

fact,  as  the  satellite  mass  comes  closer  to  that  of  the  tip  mass,  the  center  of  mass  of  the 
total  system  will  move  along  the  appendage  length.  One  example,  where  the  satellite 
mass  and  tip  mass  are  equal,  places  the  center  of  mass  equidistant  from  the  satellite  and 
tip  mass.  In  this  situation,  the  position  vector  of  the  system  center  of  mass  in  the  body 
frame  is  time  varying  as  well  as  the  system’s  moments  of  inertia.  Now,  when  a  torque 
is  generated  at  the  tip  mass,  the  controlling  body  will  experience  an  equal  torque  in  the 
opposite  direction. 

A  concern  arises  when  the  hybrid  coordinate  equations  (HCEs)  (Eq.  2.27)  are  used 
for  the  free-free  system  without  modification.  The  coupling  term  is  the  representation 
of  how  torque  between  the  controlling  body  and  the  flexible  appendage  is  transferred 
at  the  connection  point.  If  the  HCEs  are  unmodified  when  applied  to  the  free- free 
boundary  condition,  the  system  plant  demonstrates  responses  similar  to  that  of  an 
inverted  pendulum  where  one  or  more  open  loop  poles  are  placed  on  the  real  axis  and 
mirrored  in  the  left  and  right  half  parts  of  the  plane.  The  eigenvalues  for  the  system 
plant  which  includes  the  first  three  resonant  modes  are  listed  in  Table  2.5. 

A  way  to  visualize  the  impact  of  not  modifying  the  HCEs  for  the  free-free  condi¬ 
tion  is  to  look  at  the  open  loop  response  of  the  system.  If  a  gravity  gradient  spacecraft 
orientation  is  initially  roll=5°,  pitch=0°,  and  yaw=0°,  the  open  loop  response  should 
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show  the  spacecraft  oscillating  between  ±5°  if  no  energy  is  lost  from  the  system.  How¬ 
ever,  this  is  not  the  case  (see  Figure  2.10). 
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Figure  2.10:  Open  Loop  Response  of  Unmodified  HCEs 
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If  a  flexible  system  is  modeled  correctly,  it  should  contain  a  combination  of  rigid 
body  and  flexible  modes.  The  rigid  body  dynamics  will  generate  \  terms  while  the 
flexible  dynamics  produce  terms.  This  places  the  open  loop  poles  either  at 

the  origin  or  very  close  to  the  imaginary  axis  for  lightly  damped  poles.  The  eigenvalues 
listed  in  Table  2.5  indicates  the  system  model  is  incorrect. 

Consider  an  example  where  a  designer  is  attempting  to  model  an  undamped 
resonant  mode  with  a  natural  frequency  of  1  rad/sec.  The  designer  derives  the  following 
two  system  plants: 


Case  1 


V1 


1 

0 


\ 

/ 


/ 

Case2  = 

V 

which  produce  characteristic  equations  of  s2 


(2.40) 


-1  0  J 

-  1  for  Casel  and  s2  +  1  for  Case2.  The 
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Figure  2.11:  Pole-Zero  Map  for  [ 


pole  placement  for  both  cases  are  shown  in  Figure  2.11  and  Figure  2.12.  A  sign  error  in 
the  system  plant  can  change  the  desired  model  for  Case2  to  that  of  an  unstable  system 
shown  in  Casel. 

The  choice  of  mode  shapes  used  to  describe  the  deformation  of  the  system  is 
crucial.  The  cantilever  modes  apply  to  a  restrained  system  in  which  the  controlling 
body  is  assumed  to  neither  translate  nor  rotate  when  the  mode  shape  matrix  is  formed. 
A  better  assumption  is  to  no  longer  restrain  the  system  and  use  the  free-free  mode 
shapes  in  the  equations  of  motion.  This  allows  complete  freedom  of  motion  in  response 
to  impressed  moments  and  forces.  In  a  Tisserand  frame[38],  the  expressions  for  angular 
momentum  and  kinetic  energy  are  structurally  simplified  by  moving  the  axes  so  as  to 
set  the  internal  angular  momentum  always  to  zero.  The  requirement  also  makes  the 
internal  linear  momentum  zero.  This  constraint  is  accomplished  by  locating  the  origin 
of  the  frame  at  the  center  of  mass  of  the  system.  Now  that  the  frame  moves  with  the 
body,  a  floating  reference  frame,  the  measured  displacements  relative  to  this  frame  will 
be  small  where  an  inertial  frame  will  see  large  displacements  as  the  body  undergoes 


58 


Pole-Zero  Map 


1  - 

0.8 

0.6  - 

0.4  - 

■r  0.2 
< 

n 

«  o  - 

'6> 

ro 

E  -0.2  - 
-0.4 
-0.6  - 
-0.8  - 
-1 


-1 


-0.8  -0.6  -0.4  -0.2  0  0.2  0.4  0.6 

Real  Axis 


Figure  2.12:  Pole-Zero  Map  for 
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appreciable  rotation. 

The  use  of  the  floating  reference  frame  is  shown  in  Buttrill’s  study[31]  for  flexible 
aircraft  but  the  methods  are  valid  for  a  spacecraft  as  well.  The  following  assumptions 
used  in  the  study  will  also  apply  to  this  research: 

•  the  spacecraft  is  idealized  as  a  collection  of  lumped  mass  elements,  each  being 
a  finite  rigid  body,  and  each  having  an  associated  mass  and  moments  of  inertia 

•  the  elastic  restoring  force  resulting  from  displacement  of  any  mass  element  is 
linear  and  proportional  to  that  displacement 

•  the  total  rotational  displacement  of  any  lumped  mass  with  respect  to  its  unde¬ 
formed  orientation  is  small 

•  deformation  is  described  by  a  linear  sum  of  mode  shapes  multiplied  by  their 
time  dependent  participation  coefficients 

the  implications  of  the  above  assumptions  are  summarized  in  the  following  statements: 

•  each  mass  element  resides  at  a  node  of  the  structural  finite-element  model  and 
constitutes  a  lumped  resistance  to  acceleration 

•  proportional  strain  to  stress  relationships 

•  tip  deflection  is  <  10%  of  beam  length 

•  u  =  as  shown  in  Eq.  2.20 

To  satisfy  the  minimum  relative  kinetic  energy  requirement  of  a  floating  reference 
frame,  Bucken’s  constraint  relationships  must  be  considered[103]. 

Tii^ijdrrii  =  0 

^ifi<j>ijdmi  =  0  (2-41) 
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if  the  assumed  mode  shapes  are  the  eigenvectors  of  a  structural  model  in  undamped 
vibration  with  free-free  boundary  conditions,  then  motion  according  to  each  mode  shape 
j  should  satisfy  the  conditions. 

These  constraints  define  the  assumed  modes  method  integral  terms,  X\  and  A3, 
shown  in  Eq.  2.29.  When  applied  to  Eq.  2.28, 


8  =  XfR-Xj 

8  =  0 


which  reduces  Eq.  2.27  to 


16  =  Tc 

V  +  2(WnV  +  ^nV  =  0 


(2.42) 


(2.43) 


and  there  no  longer  exists  angular-vibrational  momentum  exchange. 

Kakad  provides  a  solution  which  does  not  rely  on  the  hybrid  coordinate  equation  [103]. 
This  study  discusses  the  dynamics  and  control  of  slewing  maneuvers  of  a  large  flexible 
spacecraft  named  NASA  Spacecraft  Control  Laboratory  Experiment  (SCOLE).  The 
system  was  modeled  as  a  distributed  parameter  beam  with  two  end  masses.  The  three 
dimensional  linear  variation  analysis  of  this  free-free  beam  model  is  incorporated  to¬ 
gether  with  rigid-slewing  maneuver  dynamics.  Beam  vibrations  at  the  end  of  slewing 
maneuvers  were  controlled  using  the  infinite  time  Linear  Regulator  Problem  formula¬ 
tion.  The  results  illustrated  how  slew  angle  changes  vs  time  changed  between  a  model 
using  only  rigid  body  modes  and  a  model  which  added  the  first  two  flexible  modes  to 
the  rigid  body  equations. 

Kakad’s  analysis  of  higher  uncontrolled  modes  indicated  serious  control  spillover 
due  to  coupling  among  the  modes.  This  indicated  the  residual  modes  were  excited  by 
feedback  control  which  is  designed  for  a  low-order  model,  senses  and  actuates  higher- 
order  modes,  and  renders  the  system  unstable[14]  [124].  Spillover  can  be  avoided  by 
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reducing  the  control  gain.  However,  this  often  results  in  poor  performance.  To  improve 
performance,  control  objectives  lead  to  high-gain,  high-bandwidth  control  strategies 
which  reduce  inherent  frequency  separation  between  rigid  body  and  elastic  dynamics 
and  increases  the  possibility  for  adverse  coupling  when  control  loops  are  closed [233]. 

The  SCOLE  system  modeled  by  Kakad  deals  with  slewing  maneuvers  of  an  an¬ 
tenna  attached  at  one  end  of  a  shuttle  robotic  arm.  Although  similar  techniques  were 
used  by  Kakad  in  his  analysis,  free-free  modeling  with  linear  quadratic  control,  several 
differences  between  that  work  and  the  research  presented  in  this  document  exist.  The 
on-orbit  mass  of  the  shuttle  is  approximately  100,000  kg  while  the  spacecraft  bus  mod¬ 
eled  here  is  40  kg.  Kakad’s  results  demonstrated  the  change  in  slew  response  between 
a  rigid  body  model  and  one  which  includes  the  first  two  flexible  modes.  However,  he 
doesn’t  mention  the  robustness  of  the  system,  what  the  cost/benefits  are  of  including 
or  not  including  flexible  modes,  nor  provides  an  explanation  as  to  why  he  included  the 
first  resonant  modes  and  not  higher  order  frequencies.  At  no  point  in  his  paper  does  he 
mention  the  ’’rule  of  thumb”  of  considering  the  system  as  a  rigid  body  (information  on 
crossover  frequency  in  relation  to  the  first  resonant  mode  was  not  provided). 

While  Kakad’s  study  considered  a  free-free  system,  it  was  similar  to  several  other 
studies  focused  on  flexible  dynamic  control  of  a  rotating  hub  with  an  attached  whip 
appendage  [242]  [248].  These  studies  focused  on  single  axis  analysis  to  avoid  the  adverse 
coupling  issues.  Although  the  systems  were  not  large  flexible  spacecraft,  the  mass  of 
the  hub  was  still  much  greater  than  the  mass  of  the  appendage.  Also,  the  hub  was  fixed 
in  translation  which  further  limited  the  center  of  mass  to  small  variations.  Simulation 
and  comparison  studies  done  by  Guo-Ping  show  that  even  small  tip  masses  may  affect 
dynamic  characteristics  of  the  system  significantly,  which  may  result  in  the  largening  of 
vibrating  amplitude  and  the  descending  of  vibrating  frequency  of  the  beam,  and  may 
affect  end  position  of  the  hub-beam  system  as  well[32].  While  the  efforts  are  closer  to 
the  approach  provided  within  this  research  document,  they  only  looked  at  single- input 
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single-output  systems. 

A  recent  study  which  uses  finite  element  methods  similar  those  presented  within 
this  document  was  done  by  LoBosco  in  his  modeling  of  a  large  optical  telescope  used 
to  find  terrestrial  planets  in  other  solar  systems[132].  The  optical  telescope  modeled 
was  a  large  flexible  structure  constructed  from  cantilever  beams.  A  concern  identified 
in  the  results  of  the  study  pertained  to  how  extremely  sensitive  the  optical  performance 
was  to  truss  uncertainties.  It  was  presented  the  cause  of  the  sensitivity  rests  within 
the  reaction  wheel  assembly  broadband  disturbance  and  further  study  in  this  area  was 
identified.  The  error  may  lie  in  the  fact  that  cantilever  mode  shapes  were  used  in 
the  model  instead  of  free-free  mode  shapes.  While  the  optical  telescope  has  typical 
characteristics  of  a  large  flexible  structure,  it’s  highly  precise  performance  requirements 
may  be  better  met  if  a  modeling  approach  mirroring  the  techniques  presented  within 
this  document  were  adopted. 

Given  the  small  flexible  nature  of  the  system  modeled  in  this  research  effort, 
and  the  areas  of  concern  identified  from  previous  researchers  in  the  area  of  flexible 
spacecraft  control,  the  FEM  analysis  outlined  in  Section  2.2  was  modified  for  the  free- 
free  assumption.  Material  properties  were  kept  the  same  and  the  free-free  analysis 
generated  eighteen  mode  shapes.  The  applicable  modes  and  the  frequencies  at  which 
they  occur  are  listed  in  Table  2.6  while  plots  of  the  first  three  bending  modes  and  the 
torsional  mode  are  provided  in  Figures  2.13-2.16. 

All  of  the  resonant  modes  for  the  free- free  condition  occur  at  similar  frequencies  as 
those  of  the  cantilever  modes  except  for  the  first  bending  mode.  Where  the  first  bending 
mode  of  the  cantilever  appendage  occurs  at  0.3947  Hz,  the  first  bending  mode  of  the 
free-free  appendage  is  almost  an  order  of  magnitude  higher  at  3.1518  Hz.  Therefore, 
each  assumption  will  generate  differing  mode  shape  matrices. 
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Table  2.6:  Free-Free  Resonant  Frequencies  From  FEM  Analysis 


Mode  Shape 

Frequency  (Hz) 

1st  Bending  in  Y  direction 

3.1518 

1st  Bending  in  X  direction 

3.1518 

1st  Torsional  Mode 

8.0284 

2nd  Bending  in  Y  direction 

12.588 

2nd  Bending  in  X  direction 

12.588 

3rd  Bending  in  Y  direction 

29.491 

3rd  Bending  in  X  direction 

29.491 

1st  Compression  Mode 

55.555 

4th'  Bending  in  Y  direction 

61.690 

4th  Bending  in  X  direction 

61.690 

5th  Bending  in  Y  direction 

115.70 

5th  Bending  in  X  direction 

115.70 

6th  Bending  in  Y  direction 

188.62 

6th  Bending  in  X  direction 

188.62 

7th  Bending  in  Y  direction 

279.30 

7th  Bending  in  X  direction 

279.30 

8th  Bending  in  Y  direction 

386.99 

8tfl  Bending  in  X  direction 

386.99 
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2.3  State  Space  Form  of  Nominal  System  Model 


The  FEM  analysis  output  is  a  powerful  solution  to  producing  the  equations  of 
motion  for  a  small  flexible  spacecraft  using  free-free  mode  shapes.  To  illustrate  this, 
consider  a  spacecraft  using  the  material  properties  shown  in  Section  2.2  where  the  FEM 
model  is  comprised  of  101  nodes  with  the  controlling  body  treated  as  a  lumped  mass 
located  at  node  1  and  the  tip  mass  is  a  lumped  mass  located  at  node  101.  Each  node 
is  free  to  both  translate  and  rotate  such  that  the  state  vector,  x,  takes  a  similar  form 
as  Eq.  2.25  but  is  for  101  nodes  instead  of  3.  Since  each  node  has  6  degrees  of  freedom, 
the  dimensions  of  x  is  (606  x  1). 

The  harmonic  equation  shown  in  Eq.  2.21  is  now  written  as 


[m]x  +  [k]x  =  F  (2.44) 

where  \m\  and  [k]  are  the  mass  and  stiffness  coefficients  coupling  each  node  and  F  is  the 
matrix  of  external  forces  applied  to  the  system. 

The  state  space  representation  of  Eq.  2.44  requires  a  set  of  first  order  equations. 
This  is  done  in  the  following  manner: 


let 


x  +  [m]  L[k]x  =  [m]  lF 


x  = 


=  —  [m]  L[k\x  +  [m]  lF 


(2.45) 


Zl  =  x 

Z\  =  X 

Z2  =  Zl 


Z2  =  —[m]  1  [k]  zi  +  [m]  lF 


(2.46) 


then 
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yz2  j 
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m~1F 


(2.47) 
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For  this  system  model,  F  is  the  control  torque  moments  acting  on  the  controlling 
body  located  at  node  1.  Since  there  are  three  independent  reaction  wheels  providing 
the  control  input,  F  takes  the  form 


'  03.  ' 


V 


) 


F=  I  4x3  I  «  (2.48) 

0600x3 

Attitude  measurements  for  this  system  are  the  rotational  orientations  of  the  con¬ 
trolling  body  and  is  represented  as 


y  =  Cz 


(2.49) 


where 

C  =  [03x3^3x303x60003x606]  (2.50) 

since  the  measurements  are  the  displacement  rotations  with  no  rate  measurements. 

Now,  the  state  space  representation  of  the  system  using  physical  coordinates,  x, 
can  become  overwhelming  depending  on  the  number  of  nodes  used  in  the  FEM.  For  this 
case,  the  plant  matrix,  A,  has  dimensions  (1212  x  1212).  Computer  processing  capacity 
may  be  of  little  concern  when  using  the  latest  desktop  technology,  but  such  luxury  is 
typically  not  available  to  small  satellites.  Therefore,  the  modal  approach,  shown  in 
Eq.  2.22  and  Eq.  2.23,  is  desirable  due  to  a  reduced  number  of  states  and  is  written  as 


rj  +  Df]  +  Krj  =  F 


(2.51) 


where  the  bar  denotes  a  generalized  matrix  in  modal  equations,  structural  damping  is 
included  in  the  form 

D  =  2  C,\fR  (2.52) 

and  the  modal  forces  produced  by  controller  inputs  are 


F  =  [</)}tF 


(2.53) 
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The  same  development  of  Eq.  2.47  is  applied  to  Eq.  2.51  to  generate  the  state 
space  representation  of  the  modal  equations. 

\ 


( ,  \ 

( 

Zi 

\h) 

V 

0  I 
-K  -D 


\Z2  J 


+ 


\FJ 


(2.54) 


and 


V  =  C[</>\z 


(2.55) 


Consider  what  occurs  when  this  approach  is  applied  to  cantilever  boundary  condi¬ 
tions.  The  modal  force  vector,  with  momentum  actuators  placed  only  at  the  controlling 
body,  takes  the  form 


F  =  [4>]tF 


1  N1T1MS1  N1T2MS1 
N1T1MS2 


^  NITIMSm 


NnRSMSl 


NnRZMSm  J 


0  0 
0  0 
1  0 
0  1 
0  0 
0  0 


0 

0 

0 

0 

1 

0 


N1R1MS1 

N1R2MS1 

N1R3MS1 

N1R1MS2 

N1R2MS2 

N1R3MS2 

y  0  0  0  J 


(2.56) 


^  NIRlMSm  N\R2MSm  NlR3MSm  J 


For  the  restrained  system  using  cantilever  mode  shapes,  as  used  extensively  in  the 
literature,  the  boundary  conditions  applied  to  the  first  node  constrain  both  translation 
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and  rotation.  This  means  there  is  no  rotation  of  node  1  and 


F  = 


(2.57) 


Thus,  for  the  cantilever  assumption,  the  described  method  will  generate  linearly  in¬ 
dependent  modal  equations  but  with  no  controlling  forces.  This  is  perhaps  why  the 
analytical  solutions  presented  in  the  literature,  those  originally  derived  from  Likins,  all 
assume  cantilever  modes;  which  is  a  valid  assumption  for  cases  where  the  system  center 
of  mass  is  located  at  the  center  of  mass  of  the  controlling  body  and  moves  very  slightly 
from  that  location  over  time  as  seen  with  large  space  structures  or  symmetric  spacecraft. 

Recall,  the  dimensions  of  A  from  Eq.  2.47  is  (1212  x  1212).  Now,  the  size  of 
A  is  dependent  on  the  number  of  mode  shapes  included  in  the  system  model.  There 
are  an  infinite  number  of  mode  shapes  which  can  be  generated  for  a  flexible  system. 
It  is  impossible  to  include  all  of  the  mode  shapes  since  the  accuracy  which  Patran  can 
determine  these  mode  shapes  degrades  at  higher  frequencies  and  the  mode  shape  matrix, 
[(/>],  no  longer  is  orthonormalized  if  the  number  of  modes  exceeds  the  number  of  nodes 
used  in  the  FEM.  A  control  designer  can  use  the  following  technique  to  determine  how 
many  mode  shapes  to  include  in  the  nominal  plant  model. 

The  higher  the  frequency  of  a  mode  shape,  the  smaller  the  energy  which  is  stored 
in  that  shape.  Harmonic  motion  of  the  flexible  structure  will  subject  the  controlling 
body  to  harmonic  excitation.  The  nondimensional  ratio,  ^ ,  is  a  measure  of  the  force 
transmitted  to  the  controlling  body  and  is  written  as  [147] 

A  =  [1  +  (^)2]J|G(W)|  (2.58) 

where  Q  is  the  viscous  damping  factor,  iv  is  the  excitation  frequency,  t on  is  the  natural  fre¬ 
quency  of  undamped  oscillation,  and  \G(iu>)\  is  the  magnitude  of  the  system’s  frequency 
response.  For  higher  frequency  resonant  modes,  the  force  transmitted  to  the  controlling 
body  decreases.  In  addition,  the  gain  of  a  stable  system  drops  off  at  frequencies  higher 
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Figure  2.17:  Nondimensionalized  Resonant  Transmission  Forces 


than  the  control  dynamics  of  the  closed  loop  system. 

Using  the  free-free  resonant  frequencies  listed  in  Table  2.6,  a  plot  is  created  of 
■\  versus  the  resonant  mode  number  where  each  mode  number  represents  each  mode 
shape,  with  both  bending  modes  in  the  x  and  y  direction  receiving  one  number  instead 
of  two  separate  numbers  (see  Figure  2.17). 

Resonant  modes  beyond  the  fourth  bending  mode  will  be  excluded  from  the  nom¬ 
inal  model  because  the  slope  of  the  plot  approaches  zero  beyond  that  point  which  indi¬ 
cates  the  amount  of  force  transmitted  at  higher  frequencies  are  indistinguishable  from 
one  another.  Table  2.7  lists  the  resonant  modes  generated  in  the  FEM  analysis,  which 
are  used  to  build  the  nominal  system  model.  These  modes  lead  to  a  modal  A  matrix  of 
dimensions  (32  x  32),  which  is  2.6%  the  size  of  the  plant  matrix  generated  using  physical 


coordinates. 
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Table  2.7:  Free-Free  Resonant  Frequencies  for  Nominal  Model 


Mode  Shape 

Frequency  (Hz) 

Low  Frequency  Mode 

8.3781  xl(T7 

Low  Frequency  Mode 

3.2062x  10_ti 

Low  Frequency  Mode 

3.9299x  10-e 

Low  Frequency  Mode 

1.0974x  10-5 

Low  Frequency  Mode 

1.2097xl0"5 

Low  Frequency  Mode 

1.2959x  10-5 

1st  Bending  in  Y  direction 

3.1518 

1st  Bending  in  X  direction 

3.1518 

1 st  Torsional  Mode 

8.0284 

2nd  Bending  in  Y  direction 

12.588 

2nd  Bending  in  X  direction 

12.588 

3rd  Bending  in  Y  direction 

29.491 

3rd  Bending  in  X  direction 

29.491 

1st  Compression  Mode 

55.555 

4th  Bending  in  Y  direction 

61.690 

4th  Bending  in  X  direction 

61.690 

Chapter  3 


Attitude  Control  of  Small  Flexible  Structures 


3.1  Control  Bandwidth 

The  normal  mode  solution  of  the  FEM  model  generates  modal  information  for  all 
resonant  modes  of  the  flexible  structure.  While  all  of  this  information  is  used  in  creating 
the  nominal  model  of  the  system,  the  design  model  may  have  fewer  resonant  modes 
included.  The  stability  and  robustness  of  an  optimal  controller  may  be  used  by  designers 
to  determine  which  modes  to  include  (see  Section  3.3).  Knowledge  of  the  proximity  of 
resonant  modes  near  the  controller  bandwidth  is  also  an  important  consideration.  The 
mode  frequencies  near  the  bandwidth  of  the  satellite’s  attitude  control  system  will  have 
a  larger  impact  on  control  authority  than  the  higher  frequency  modes.  According  to 
some  designers,  the  lowest  natural  frequencies  of  flexible  components  should  be  at  least 
an  order  of  magnitude  greater  than  the  rigid-body  frequencies  before  flexibility  can  be 
neglected[3].  This  simplification  of  design  is  only  valid  in  models  based  on  assumed 
cantilever  boundary  conditions  and  has  not  been  validated  for  system’s  using  free-free 
boundary  conditions. 

The  controller  bandwidth  defines  the  frequency  where  the  control  authority  begins 
to  diminish.  Attitude  control  and  disturbance  rejection  are  effective  from  DC  up  to 
the  bandwidth [235]  and  high  accuracy  implies  high  position  gain  and  high  bandwidth. 
However,  increasing  bandwidth  may  cause  bending  resonances  to  affect  control  system 
performance  by  producing  large  output  errors  due  to  measurement  noise. 
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Control  designers  can  estimate  controller  closed  loop  bandwidth  with  the  following 
equation: 


(jJBW  = 


(3.1) 


where  I\p  is  the  position  gain  of  the  controller  and  I  is  the  spacecraft’s  moment  of 
inertia.  The  position  gain  can  be  estimated  and  must  be  high  enough  to  provide  the 
required  attitude  control  pointing  accuracy  in  the  presence  of  disturbance  torques.  The 
units  of  Kp  can  either  be  or  and  is  determined  by 


K  >  — 

p~  9f 


(3.2) 


where  To  is  the  peak  disturbance  torque  and  6e  is  allowable  attitude  error. 

Four  environmental  disturbance  torques  are  commonly  considered  when  determin¬ 
ing  the  peak  disturbance  torque.  These  were  the  gravity  gradient,  Tg,  solar  radiation 
pressure,  Tsp,  magnetic  field,  Tm,  and  aerodynamic  torques,  Ta.  The  equations  used  to 
calculate  these  disturbance  torques  are: 

^  =  ^|4-4lsin(2^)  (3.3) 


where  Tg  is  the  max  gravity  torque,  /i  is  the  Earth’s  gravitational  constant  (3.986  x 
1014^),  R  is  the  orbit  radius  (m),  9  is  the  maximum  deviation  of  the  Z-axis  from  local 
vertical  in  radians,  and  Iz  and  Iy  are  the  moments  of  inertia  about  z  and  y  (or  x,  if 
smaller)  axes  in  kgm2.  The  solar  radiation  pressure  can  be  represented  as 


Tsp  =  F(cps  -  cg)  (3.4) 

where 

F  =  —  As(\  +  q)  cosi  (3-5) 

c 

and  Fs  is  the  solar  constant,  1,367^-,  c  is  the  speed  of  light,  3  x  108^,  As  is  the  surface 
area  (m2),  cps  is  the  location  of  the  center  of  solar  pressure  (to)  ,  cg  is  the  center  of 


74 

Table  3.1:  Typical  Disturbance  Torques  for  LEO  Small  Satellites 


Disturbance  Torque 

Worst  Case  Magnitude 

Magnetic  Torque  (D  =  1.0) 

5.1468X  10”5 

Gravity  Gradient  Torque  (Iy  —  Iz  =  139.3) 

4.3686X  10”5 

Drag  Torque  (cpa  —  cg  =  3 cm) 

3.489  xl0-y 

Solar  Pressure  Torque  ( cps  —  cq  =  3 cm) 

0.02205x10'° 

gravity  (m),  q  is  the  reflectance  factor,  and  i  is  the  angle  of  incidence  of  the  Sun.  The 
magnetic  held  torque  on  the  spacecraft  is 

Tm  =  DB  (3.6) 

where  D  is  the  residual  dipole  of  the  vehicle  {Am2)  and  B  is  the  Earth’s  magnetic  held 
in  tesla.  B  can  be  approximated  as  for  a  polar  orbit  to  half  that  at  the  equator.  M 
is  the  magnetic  moment  of  the  Earth,  7.96  x  10 15tesla  ■  m3,  and  R  is  the  radius  from 
the  dipole  center  to  the  spacecraft  in  m. 

Finally,  the  aerodynamic  drag  on  the  satellite  can  be  written  as 

Ta  =  F(cpa  -  cg)  (3.7) 

where  F  =  OApCjAV2 ,  Cd  is  the  drag  coefficient,  p  is  atmospheric  density  (-^-),  A  is 
the  surface  area  (m),  V  is  the  spacecraft  velocity  (^),  cpa  is  the  center  of  aerodynamic 
pressure  (to),  and  cg  is  the  center  of  gravity  (m). 

For  a  typical  Low  Earth  Orbit  (LEO)  small  satellite  with  a  deployed  appendage, 
the  worst  case  magnitude  for  these  environmental  disturbance  torques  are  calculated 
and  listed  in  Table  3.1  [231]: 

The  peak  disturbance  torque,  To.  used  in  Eq.  3.2  would  be  if  all  of  the  envi¬ 
ronmental  disturbance  torques  were  acting  at  the  same  time  in  resonance.  Assuming 
To  =  1  x  10“4  and  a  attitude  pointing  accuracy  requirement,  6e  =  0.1°,  the  controller 
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closed  loop  bandwidth  would  be  c obw  =  0.45Hz.  As  the  altitude  of  small  satellites 
decrease,  an  increase  in  atmospheric  drag  occurs.  If  the  same  satellite  is  placed  in 
an  altitude  comparable  to  that  of  the  space  shuttle,  the  atmospheric  density  becomes 
approximately  two  orders  of  magnitude  larger.  This  will  cause  the  drag  disturbance 
torque  to  be  approximately  2x  10~4_/Vm  and  the  control  bandwidth  to  approximately 
equal  0.7Hz.  The  flexible  nature  of  the  satellite  system  further  complicates  the  deter¬ 
mination  of  a  specific  value  for  control  bandwidth  as  the  center  of  mass  of  the  system 
relative  to  the  center  of  pressure  is  continuously  changing. 

3.2  LQG/LTR 

Linear  Quadratic  Regulator  (LQR)  control  leads  to  linear  control  laws  that  are 
easy  to  implement  and  analyze.  The  system  being  controlled  is  assumed  to  be  at 
equilibrium  and  it  is  desired  to  maintain  the  equilibrium  despite  disturbances. 


3.2.1  Linear  Quadratic  Regulator  (LQR) 

If  a  system  is  represented  as 

x  =  Ax  +  Bu 

y  =  Cx  (3-8) 


a  cost  function  can  be  defined  as 


(. x'Qx  +  vl Ru)dt 


(3.9) 


where  J  is  minimized  with  respect  to  the  control  input  u(t). 

J  represents  the  weighted  sum  of  energy  of  the  state  and  control.  Q  and  R  are 
weighting  matrices,  or  design  parameters,  where  the  state-cost  matrix,  Q,  weights  the 
states  while  the  performance  index  matrix,  R,  weights  the  control  effort.  If  Q  is  increased 
while  R  remains  constant,  the  settling  time  will  be  reduced  as  the  states  approach  zero 
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at  a  faster  rate.  This  means  that  more  importance  is  being  placed  on  keeping  the  states 
small  at  the  expense  of  increased  control  effort.  If  R  is  very  large  relative  to  Q,  the 
control  energy  is  penalized  very  heavily.  This  physically  translates  to  smaller  motors, 
actuators,  and  amplified  gains  needed  to  implement  the  control  law. 

Pontriagin’s  minimum  principle  is  used  to  solve  the  optimal  control  problem. 
First,  from  the  Hamiltonian 

H(x ,  A,  t)  =  -(x'Qx  +  u'Ru )  +  A  '{Ax  +  Bu)  (3.10) 

the  minimum  principle  states  that  the  optimal  control  and  state  trajectories  must  satisfy 
the  following  three  equations: 


x 

-A 

dH 

du 


dH 

~d\ 

8H 

dx 

o 


(3.11) 


where  x(0)  =  xo  are  the  state  equations  and  A (T)  =  0  are  the  costate  or  adjoint 
equations. 

When  Eqs.  3.8-3.11  are  applied,  the  equations  become 


x  =  Ax  +  Bu 
—A  =  Qx  +  A'  A 

u*  =  -R^B'X  (3.12) 


where  u*  is  the  optimal  control  and  R  has  to  be  positive  definite  for  its  inverse  to 
exist.  The  above  coupled  linear  differential  equations  form  a  two  point  boundary  value 
problem  (TPBVP),  which,  because  of  mixed  boundary  conditions,  is  difficult  to  solve 
numerically.  Substituting  the  optimal  control  into  the  state  equation  produces 
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u(t)  =  —K(t)x(t)  (3.18) 

where  K(t )  =  R~1B' P(t). 

For  the  infinite  time  LQR  problem,  we  let  t  approach  infinity.  Of  course,  now  one 
runs  into  the  question  of  the  convergence  of  the  cost  function  and  the  existence  of  the 
optimal  controller.  Even  if  the  optimal  control  exists,  it  does  not  necessarily  result  in  a 
stable  closed  loop  system.  It  turns  out  that  under  mild  conditions,  P(t)  approaches  a 
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constant  matrix  P  — ? >  0)  and  the  positive  definite  solution  of  the  algebraic  Riccati 

equation  results  in  an  asymptotically  stable  closed  loop  system. 

Q  +  A'P  +  PA-  PBRt'B'P  =  0 

u  =  —  Kx 

K  =  R~lB'P  (3.19) 

The  exact  conditions  for  the  above  to  hold  are  the  following.  The  pair  (A,B) 

are  stabilizable,  R  >  0,  and  Q  can  be  factored  as  Q  =  C'qCq,  where  Cq  is  any  ma¬ 

trix  such  that  {Cq, A)  is  detectable.  Stabilizable  refers  to  the  uncontrollable  portions 
being  asymptotically  stable  while  detectable  refers  to  the  unobservable  portions  being 
asymptotically  stable.  These  conditions  are  necessary  and  sufficient  for  existence  and 
uniqueness  of  an  optimal  controller  that  will  asymptotically  stabilize  the  system. 

If  one  is  interested  in  controlling  a  subset  of  the  states,  the  system  outputs  for 
this  example,  the  cost  function,  J,  can  be  written  as 

1  r°° 

J  =  -  /  {y'Qy  +  u'Ru)dt 

2  Jo 

1  r°° 

=  -  /  ( x'C'QCx  +  u Ru)dt  (3.20) 

2  Jo 

and  the  Hamiltonian  takes  the  form 

'  A  -BR^B'  ^ 

H  =  (3.21) 

v  -C'QC  -A' 

An  interesting  thing  to  note  is  the  2n  eigenvalues  of  H  are  symmetric  about  both 
the  imaginary  axis  and  the  real  axis.  Thus,  the  adjoined  system  has  n  stable  roots  and 
n  unstable  roots;  half  associated  with  x  and  the  other  half  associated  with  A.  For  the 
cost  function,  J,  to  remain  finite,  the  n  stable  eigenvalues  of  H  must  be  the  closed  loop 
poles  of  the  system.  When  we  optimize  a  controllable  linear  system  using  a  quadratic 
cost,  we  will  always  generate  a  stable  closed-loop  system. 
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3.2.2  Optimal  Estimation  (Kalman  Filter) 

The  LQR  solution  is  basically  a  state  feedback  type  of  control  meaning  it  requires 
that  all  states  be  available  for  feedback.  This  is  usually  unreasonable  and  some  form  of 
state  estimation  is  necessary.  The  combination  of  the  state  feedback  and  observer  will 
always  result  in  a  stable  closed-loop  system.  Controller  performance  can  be  optimized 
according  to  some  quadratic  cost  function  (LQR).  Observer  design  can  also  be  done  in 
an  optimal  manner  provided  the  problem  is  formulated  in  a  probabilistic  (or  stochastic) 
sense.  The  formulation  of  the  state  estimation  problem  is  as  follows: 

x  =  Ax  +  Bu  +  co 

y  =  Cx  +  v  (3.22) 

where  lo  represents  random  noise  disturbance  input  (process  noise)  and  v  represents 
random  measurement  (sensor)  noise. 

It  is  assumed  that  both  noise  processes  are  unbiased  white  Gaussian  zero-mean 
stationary  processes  with  known  covariances  given  below. 

E{uj(t)}  =  0 
E{v(t)}  =  0 

E{u>(t)u(t  +  t)'}  =  Qo5{t-r) 

E{v(t)u(t  +  t)1}  =  Ro5(t  —  r) 

E{uj(t)v(t  +  t)'}  =  0  (3.23) 

The  state-space  solution  to  this  problem  was  first  introduced  by  R.  E.  Kalman 
and  R.  S.  Bucy.  It  obtains  an  estimate  of  x(t)  based  on  noise-corrupted  measurements 
such  that  the  variance  of  the  error  is  minimized.  The  optimal  estimator  is  given  by 


x  =  Ax  +  Bu  +  L(y  —  Cx) 


(3.24) 
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where  x  is  the  estimate  of  x.  The  observer  gain  is  computed  from 

L  =  TC'R^1  (3.25) 

and  X  is  found  as  the  positive  semi-definite  solution  of 

TX  +  XT'  +  Q0  -  TC'Rq lCT  =  0  (3.26) 

Note  that  the  equations  for  the  filter  gain,  L,  and  X  are  very  similar  to  the 
equations  for  the  LQR  solution.  In  particular,  the  equation  for  X  is  an  algebraic  Riccati 
equation  and  is  the  estimation  error  covariance  and  the  trace  of  X  indicates  how  well 
the  filter  is  performing.  The  Q o  and  Rq  matrices  represent  the  intensity  of  the  process 
and  sensor  noise  inputs  and  are  the  only  parameters  that  are  to  be  provided  by  the  user. 

3.2.3  Linear  Quadratic  Gaussian  (LQG) 

The  LQG  combines  the  LQR  with  the  Kalman  filter 

x  =  Ax  +  Bu  +  uj 
y  =  Cx  +  v 

with  controller 

u  =  —  Kx 

K  =  R~lB’P 

0  =  A'P  +  PA  +  Q  -  PBR~lB'P  (3.28) 

and  observer 

x  =  Ax  +  Bu  +  L{y  —  Cx) 

L  =  XC'Rp1 


.  Consider  the  plant  equations 


(3.27) 


0 


AT,  +  XT'  +  Q0  -  TC'R^CT 


(3.29) 
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The  LQR  results  in  an  asymptotically  stable  closed  loop  system.  In  addition,  the 
controller  minimizes  the  average  of  the  LQR  cost  function  (weighted  variance  of  the 
state  and  input) 

J  =  -  /  (. x'Qx  +  u' Ru)dt  (3.30) 

2  Jo 

producing  an  optimal  solution. 

The  transfer  function  of  the  LQG  compensator  is  similar  to  the  observer-based 
compensator,  and  is  given  by 

H(s)  =  K(sI-A  +  BK  +  LC)~lL  (3.31) 

An  important  thing  to  note  is  that  LQG  has  no  guaranteed  stability  margins  like 
those  produce  by  LQR.  In  fact,  its  margins  can  be  dangerously  low.  By  changing  the 
design  parameters  Q,  R,  and  the  noise  intensities,  it  is  observed  that  some  parameters 
can  have  drastic  effects  on  the  system  properties. 

3.2.4  Loop  Transfer  Recovery  (LTR) 

The  major  problem  with  the  LQG  solution  is  its  lack  of  robustness.  The  loop 
transfer  recovery  (LQG/LTR)  technique  maintains  the  LQG  machinery  but  modifies  the 
design  procedure  to  address  some  of  the  short  comings  of  the  original  LQG  approach. 
The  open  loop  transfer  function  of  the  LQR  is  given  by 

L(s )  =  K$(s)B  (3.32) 

where  <5(s)  =  (si  —  A)-1. 

The  open  loop  transfer  function  for  LQG  is  likewise  given  by 

L(s) lqg  =  K(sl  -  A  +  BK  +  LC)~1LC^(s)B  (3.33) 

Under  the  following  two  conditions 

(1)  Gp,  the  system  plant,  is  minimum-phase  (i.e.  it  has  no  zero  in  the  RHP) 
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(2)  Ro  =  1  and  Qo  =  q2BB' 
it  can  be  shown  that 

Jim  L(s)Lqg  =  L(s)  (3.34) 

This  suggests  the  following  procedure  for  design.  Choose  the  LQR  parameters  (Q 
and  R)  such  that  the  LQR  loop  transfer  function,  L(s),  also  called  the  target  feedback 
loop  (TFL),  has  desirable  time  and/or  frequency  domain  properties.  Design  an  observer 
with  parameters  specified  in  (2)  above.  Increase  the  tuning  parameter,  q,  until  the 
resulting  loop  transfer  function  is  as  close  as  possible  to  the  TFL  while  considering 
Bode  gain,  closed  loop  bandwidth,  and  control  effort  limitations.  Because  the  loop 
transfer  function  of  LQG  approaches  that  of  LQR,  it  will  asymptotically  recover  its 
properties. 

Loop  Shaping  Step 

(1)  Determine  the  controlled  variables  (which  may  or  may  not  be  the  same  as  the 
measured  variables)  and  set  Q  =  C'C  or  Q  =  C'qCq 

(2)  Convert  the  design  specifications  into  a  desired  TFL 

(3)  Vary  the  parameter,  q,  until  the  resulting  loop  transfer  function  is  similar  to 
the  TFL.  One  may  use  the  root  square  locus  (RSL)  approach  here  for  SISO 
systems  [2 12], 

To  accomplish  the  recovery  step,  select  a  scalar,  q,  and  solve  the  filter  Riccati 
equation 

AS  +  S A'  +  q2BB'  -  SC'CS  =  0  (3.35) 

and  set  L  =  SCL 

Increase  q  until  the  resulting  loop  transfer  function  is  close  to  the  TFL.  The  higher 
the  value  of  q,  the  closer  the  LQG  system  comes  to  the  LQR  performance.  It  should  be 
noted  that  the  value  of  q  should  not  be  increased  indefinitely  because  this  may  lead  to 
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unreasonably  large  values  for  the  controller  gain.  Also,  because  LQR  has  -20dB  slope 
at  high  frequencies,  large  values  of  q  will  also  recover  this  slow  roll-off  rate.  Smaller 
values  for  q  will  tend  to  trade  off  lower  stability  margins  with  higher  roll-off  rates  at 
higher  frequencies. 

3.3  Reduced-Order  Controllers 

The  basic  problem  in  controlling  a  flexible  structure  is  the  presence  of  a  large 
number  of  lightly  damped  structural  modes.  Theoretically,  there  exist  an  infinite  num¬ 
ber  of  resonant  modes.  The  higher  the  modal  frequency,  the  smaller  the  energy  which  is 
contained  within  the  mode.  It  isn’t  feasible  for  a  numerical  analysis  to  be  conducted  on 
a  nominal  system  plant  model  which  contains  one  thousand,  or  even  one  hundred,  reso¬ 
nant  modes.  The  energy  of  the  higher  frequency  modes  will  have  minimal  impact  on  the 
overall  dynamic  response  and  practical  limitations  necessitate  the  use  of  reduced-order 
controllers.  The  initial  analysis  runs  with  Patran,  shown  in  Section  2.3,  generated  16 
resonant  modes. 

3.3.1  Controller  Robustness 

The  uncontrolled  modes,  as  well  as  the  error  in  the  knowledge  of  the  controlled 
modes,  represent  uncertainty.  Since  the  number  of  structural  modes  is  usually  large  and 
finite  element  modeling  accuracy  typically  decreases  with  increasing  modal  frequency, 
the  design  model  should  contain  the  first  few  resonant  modes.  The  remaining  structural 
modes  then  constitute  the  plant  uncertainty. 

For  single-input  single-output  systems,  the  relationship  of  stability,  sensitivity  re¬ 
duction,  disturbance  attenuation  and  rejection  to  the  return  difference  of  the  closed  loop 
system  has  been  understood  for  quite  some  time[22]  [87].  Attempts  at  extending  classi¬ 
cal  design  procedures  to  multivariable  systems  concentrated  on  sensitivity  reduction [52] 
and  examining  scalar  quantities  associated  with  the  return  difference  matrix,  which 
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include  inverse  Nyquist  procedures  for  diagonally  dominant  systems  [179]  and  charac¬ 
teristic  loci  plots[137]  [136].  However,  the  designs  may  possess  undetected  sensitivities 
to  simultaneous  perturbations  in  the  elements  of  the  return  difference  matrix[57].  In¬ 
terpretations  of  quantities  such  as  gain  and  phase  margins  in  each  input  channel  which 
lead  to  invaluable  design  insights  in  the  single-input  single-output  case  do  not  extend 
to  the  multivariable  generalizations  [53]. 

A  better  understanding  of  stability  margins  for  multivariable  systems  is  found  in 
the  investigation  of  the  norms  of  the  inverse  Nyquist  matrix.  The  measure  of  stability 
is  taken  to  be  the  magnitude  of  the  smallest  plant  perturbation  which  causes  instability 
and  is  given  by  an  appropriate  matrix  norm[57].  The  most  commonly  used  norm  is  the 
2-norm.  The  2-norms  of  a  matrix  and  its  inverse  are  the  largest  singular  value  and  the 
inverse  of  the  smallest  singular  value  of  the  matrix[63]  [118].  Singular  value  analysis 
is  popular  because  the  interpretation  of  the  smallest  singular  value  of  a  matrix  is  the 
distance  between  the  matrix  and  the  nearest  singular  matrix.  Since  this  is  precisely  the 
concept  needed  to  determine  the  nearness  of  a  stable  transfer  function  to  an  unstable 
one,  its  use  as  a  measure  of  stability  robustness  is  natural. 

Understanding  how  stability  margins  are  represented  in  multivariable  systems  is 
important.  However,  one  must  also  attempt  to  reduce  the  sensitivity  of  the  closed- 
loop  system  to  plant  perturbations  as  compared  to  the  sensitivity  to  the  open-loop 
system.  A  sufficient  condition  for  the  return  difference  matrix  to  satisfy  to  achieve  the 
desired  sensitivity  reduction  has  been  developed  in  the  form  of  a  positive  definiteness 
condition[110]  [169]  [51]  [217].  This  condition  can  be  expressed  as  a  condition  on  the 
smallest  singular  value  of  the  return  difference.  Comparison  sensitivity  and  maintenance 
of  stability  both  seek  to  retain  a  qualitative  system  property  of  sensitivity  reduction 
versus  stability  under  errors  in  the  plant  model. 
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Both  the  sensitivity  function, 


S  = 


1 


1  +  GpGc 

and  the  complementary  sensitivity  function, 

GnGr 


(3.36) 


T  = 


J  C 


(3.37) 


1  +  GPGC 

of  the  closed-loop  system  play  an  important  role  in  the  stability  robustness  of  the  system. 
Stability  requirements  bound  the  outputs  for  all  bounded  disturbances  and  bounded  ref¬ 
erence  inputs.  Doyle  and  Stein [58]  applied  Nyquist  encirclement  counts  for  the  function 


det(I  +  GPGC 


(3.38) 


to  multivariable  systems  and  placed  magnitude  constraints  to  determine  the  robustness 
of  the  system  in  the  face  of  uncertainties. 

In  classical  SISO  problems,  gain  and  phase  margins  are  used  to  characterize  tol¬ 
erable  uncertainty.  This  characterization  is  generalized  to  MIMO  problems  as 


G(s)  +  AG(s)  —  [I  +  T(s)]G,(s) 
where  L(s)  is  an  arbitrary  stable  transfer  matrix  with 

a[L{ju)\  <  m(u) 


(3.39) 


(3.40) 


and  covers  simultaneous  gain,  phase,  and  direction  errors  which  are  unknown  but 
bounded  in  size.  The  bound  m(u)  indicates  the  maximum  normalized  magnitude  which 
the  model  error  can  attain.  Stability  is  maintained  in  the  presence  of  all  possible  un¬ 
certainties,  Eq.  3.39  and  Eq.  3.40,  if  and  only  if  the  complementary  sensitivity  function 
satisfies 

<  — 7-r  (3-41) 

m(uj) 

for  all  u).  Then, 


det{I  +  GpGc  +  LGpGc)  =  det(I  +  Lr)det(5"1)  >  0 


(3.42) 
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Multiplicative  uncertainty 


Figure  3.1:  Defining  Uncertainty  in  a  Model 


and  the  uncertainty  can’t  change  the  number  of  encirclements  in  the  Multivariable 
Nyquist  Criterion.  The  system  is  stable  if  the  unperturbed  system  is  stable. 

The  loop  gain  condition  is  now  written  as 

a(L)a(T)  <  1  (3.43) 

to  arrive  at  an  upper  bound  on  the  magnitude  of  the  uncertainty  multiplier,  L,  to 
guarantee  closed  loop  stability  under  perturbations.  This  condition  is  also  necessary  if 
all  possible  multivariable  perturbations  are  allowed,  since  there  would  exist  an  L  with 

^(L)  =  W)  (3'44) 

which  brings  the  determinate  above  zero  and  changes  the  encirclement  count[121]. 

Plant  uncertainty,  errors  in  the  plant  model,  can  be  represented  as  either  multi¬ 
plicative  or  additive  uncertainty  (see  Figure  3.1).  The  multiplicative  uncertainty  form 
is  preferred  in  the  literature  on  robustness  because  the  compensated  transfer  function 
has  the  same  uncertainty  representation  as  the  nominal  model.  However,  since  flexible 
structure  models  exhibit  naturally  the  additive  uncertainty  form  of  the  transfer  function 
matrix,  this  will  be  used. 

For  the  case  of  multiplicative  uncertainty,  the  closed  loop  system  is  stable  if 


(j[Lp(ju)  -  1]  <  a  [I  +  (Gp(juj)Gc(ju>))  x] 


(3.45) 
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where  Gp(s )  and  Gc(s)  are  the  design  model  plant  and  compensator  transfer  matrices, 
and  a  and  a  denote  the  largest  and  smallest  singular  values.  At  high  frequencies, 
assuming  \\Lp(juj)\\  3>  1  and  \\Gp(juj)Gc(juj) ||  <C  1,  the  above  stability  condition  yields 

~<GpGc)  <  — (3.46) 

a\ Lp) 

The  uncertainty,  or  robustness,  barrier  is  an  upper  bound  lm{u)  on  a(Lp).  The 
system  is  stable  in  the  presence  of  unstructured  uncertainties  if  a[GpGc]  <  at 

high  frequencies. 

When  the  additive  uncertainty  formulation  is  used,  a  sufficient  condition  for  sta¬ 
bility  robustness  is  given  by  [53] 

g(/7;°fe)  >  ff(AG)  (3.47) 

vyGc) 

a  detailed  derivation  of  the  stability  robustness  conditions  with  associated  theorems  and 
proofs  can  be  found  in  [53]. 

At  high  frequencies,  assuming  ||GPGC||  <C  1,  Eq.  3.47  approximately  yields 


d(Gc)  < 


1 

u(AG) 


(3.48) 


that  is,  the  compensator  must  roll  off  sufficiently  rapidly  at  high  frequencies  to  remain 
robust  in  the  face  of  unmodelled/uncertain  high  frequency  structural  modes  and  noise. 
The  main  objective  of  the  LQG/LTR  approach  is  to  first  design  a  full  state  compensator 
which  has  the  behavior  of  the  desired  loop  transfer  matrix  (i.e.  the  loop  gain  GpGc). 
Therefore,  any  loop  shaping  should  involve  the  product  GpGc  rather  than  Gc  alone  as 
in  Eq.  3.47  and  Eq.  3.48.  Assuming  Gp  is  a  square  matrix 


Gc  =  G~\GPGC) 
a(Gc)  <  a(G~l)a(GpGc) 
°{GC)  <  a-\Gp)a(GpGc) 


(3.49) 
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Substituting  Eq.  3.49  into  Eq.  3.47,  the  following  sufficient  condition  for  stability 
robustness  is  obtained: 

a(I  +  GpGc)a(Gp) 


>  a(AG) 


(3.50) 


HGPGC ) 

The  performance  of  the  closed  loop  system  depends  on  the  low  frequency  gain 
and  the  crossover  frequency  of  the  loop  transfer  matrix  GpGc\  that  is  the  behavior  of 
(r[GpGc\.  Larger  low  frequency  gain  and  crossover  frequency  indicate  better  tracking 
performance  (also  shown  in  Section  3.1).  Thus,  (i[GpGc ]  should  lie  above  the  perfor¬ 
mance  specification.  The  other  requirement  is  the  stability  robustness  in  the  presence  of 
model  uncertainties.  If  the  multiplicative  uncertainty  formulation  is  used,  the  a[GpGc ] 
plot  should  pass  under  the  robustness  barrier  cr”1  [Lp\  at  high  frequencies.  If  the  additive 
formulation  is  used,  the  robustness  condition  of  Eq.  3.50  should  be  satisfied. 


3.3.2  Pole-Zero  Cancelation 

Looking  at  the  robustness  of  the  controller  in  the  presence  of  unmodeled  higher 
frequency  modes  most  often  leads  to  lower  order  controllers  than  those  that  are  required 
to  control  the  nominal  plant.  This  is  a  result  of  the  design  model  being  lower  in  order  due 
to  the  inclusion  of  the  lower  frequency  modes  while  high  frequency  modes  are  lumped 
into  the  uncertainty  of  the  system  dynamics.  It  may  be  possible  to  further  reduce  the 
order  of  the  controller  when  pole-zero  cancelation  is  taken  into  consideration. 

Consider  a  design  model  where  resonant  modes  (shown  in  Table  2.7)  of  frequencies 
less  than  20  Hz  are  required  to  satisfy  stability  robustness.  This  will  reduce  the  order 
of  the  controller  from  32  down  to  22.  Since  the  small  flexible  spacecraft  is  a  multi-input 
multi-output  (MIMO)  system  using  three  controller  inputs  and  measuring  three  angular 
displacements,  the  controller  is  comprised  of  nine  transfer  functions  all  of  the  order  22. 
This  high  order  of  the  controller  may  be  unnecessary  if  some  of  the  controller  poles  are 
in  close  proximity  to  controller  zeros. 

There  is  no  magic  number  to  identify  if  the  controller  poles  and  zeros  are  within 
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close  proximity  of  each  other  or  not.  Instead,  a  useful  technique  in  further  reducing 
the  order  of  the  controller  is  to  evaluate  the  effect  on  dynamic  response  if  some  poles 
are  allowed  to  cancel  out  with  zeros  in  close  proximity.  The  Matlab  command  minreal, 
minimal  realization,  is  a  useful  tool  in  accomplishing  pole-zero  cancelation.  It  is  a 
straightforward  search  through  the  poles  and  zeros  looking  for  matches  that  are  within 
a  specified  tolerance [144],  The  default  tolerance  is  y/ eps ,  where  eps  is  machine  precision. 
The  tolerance  value  can  be  increased  to  force  additional  cancelations  as  long  as  the  Bode 
plots  don’t  deviate  significantly  from  the  unreduced  controller.  If  there  is  a  difference  in 
the  Bode  plots,  a  comparison  between  the  dynamic  response  of  the  reduced  controller 
and  the  controller  prior  to  pole-zero  cancelation  will  indicate  whether  the  tolerance  was 
set  too  high.  If  the  dynamic  response  shows  little  change,  this  may  indicate  that  while 
a  certain  pole-zero  cancelation  may  generate  differing  Bode  plots,  the  contribution  of 
that  particular  pole-zero  pair  is  not  significant  enough  to  alter  the  system  response 
drastically.  Use  of  this  technique  is  demonstrated  in  Section  3.4. 

3.4  Base-Line  Response 

A  Matlab  file  was  created  to  accomplish  this  research.  This  file  reads  in  the  FEM 
data  and  generates  the  mode  shape  matrix  for  both  the  nominal  and  design  model,  as 
outlined  in  Section  2.3,  as  well  as  reading  in  the  natural  frequencies  for  each  mode.  In 
addition  to  the  system  parameters  presented  in  Section  2.2.1,  a  damping  ratio  of  0.01  for 
the  EMC  material  is  assumed.  This  assumption  is  common  for  the  structural  material 
industry  when  experimental  data  is  not  available[41]  [206]. 

A  key  design  consideration  is  the  controller  effort.  Since  the  flexible  system  is 
a  gravity  gradient  appendage  for  a  small  satellite,  it  will  take  more  control  effort  to 
correct  pitch  and  roll  displacements  than  yaw  displacements.  In  this  case,  as  well  as 
for  most  flexible  appendages,  greater  control  effort  is  required  to  produce  the  same 
dynamic  response  one  would  experience  from  a  three-axis  stabilized  cube.  The  small 
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satellite  limitations  of  mass,  volume,  and  available  power  impact  the  torque  produced 
from  on  board  reaction  wheels.  The  typical  maximum  torque  produced  by  reaction 
wheels  on  small  satellites  is  in  the  0.02Nm  to  0.3Nm  range[62][lll][157][251].  It  is 
possible  that  in  the  near  future  more  power  to  be  made  available  to  on-board  control 
actuators  through  improved  power  designs  and  increased  solar  array  efficiency.  For  this 
study,  the  maximum  allowable  control  effort  for  the  analysis  of  the  system  will  be  3Nm. 

Additionally,  the  mission  of  this  small  flexible  spacecraft  is  to  meet  the  attitude 
requirements  of  the  scientific  payloads  presented  in  Appendix  D.  The  attitude  control 
system  needs  to  keep  these  payloads  pointed  ±5°  in  the  ram,  or  velocity,  direction.  This 
requirement  applies  to  both  yaw  and  pitch  while  no  hard  requirement  is  placed  on  the 
roll  of  the  controlling  body.  Discussions  with  payload  scientists  identified  a  performance 
metric  which  concludes  that  while  data  is  usable  within  5°  off  nominal,  the  accuracy  of 
the  data  degrades  the  further  away  from  0°  the  payloads  are  pointed [208]. 

A  controller  is  designed,  using  the  LQG/LTR  techniques  provided  in  Section  3.2, 
to  optimally  reach  nominal  pointing  requirements  once  pitch/yaw  are  5°  off  nominal 
while  operating  within  the  control  effort  limit  of  3Nm.  This  initial  displacement  is  a 
result  of  the  reaction  wheels  periodically  reaching  saturation  and  momentum  dumping 
is  used  in  conjunction  with  despinning  the  wheels.  The  scientists  are  concerned  with 
how  much  time  is  required  for  the  payloads  to  go  from  5°  to  0°  and  reach  a  steady  state 
such  that  vibrations  induced  by  the  flexible  appendage  have  minimal  impact  on  system 
dynamic  response. 

Initial  values  used  for  the  design  parameters  Q  and  R  for  the  LQR  calculation, 
Q o  and  Rq  for  the  Kalman  filter,  and  q  for  the  LTR  calculation  are  shown  in  Table  3.2. 
The  design  model  for  the  initial  run  was  set  to  the  nominal  model  since  the  first  step 
is  to  determine  if  the  required  control  effort,  resulting  from  the  design  parameter  gains, 
falls  within  the  limit  of  3Nm. 

A  simulation  was  run  using  Simulink  to  plot  the  system  response  to  the  LQG/LTR 
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Table  3.2:  Initial  Values  for  Design  Parameters 


Design  Parameter 

Value 

Q 

Inxn 

R 

CO 

X 

CO 

Qo 

variance  of  1 

Ro 

variance  of  0.01 

q 

[1  le2  le4  le6] 

Max  Control  Effort 

3Nm 

controller  (see  Figure  3.2).  While  physical  coordinates  are  easier  to  visualize,  the  con¬ 
troller  design  is  based  on  modal  coordinates  because  of  the  greatly  reduced  number  of 
states.  Typically,  the  desired  attitude  is  set  to  zero  while  the  undeformed  system  has 
some  initial  orientation.  Initial  orientation  is  commonly  thought  of  in  physical  coordi¬ 
nates  with  some  initial  angular  displacement.  Keep  in  mind  the  form  of  the  physical 
coordinates. 


/ 


x  = 


N1T1  ^ 

N1T2 

N1T3 

N1R1 

N1R2 

N1R3 

N2T1 


(3.51) 


y  NnR3  J 

If  one  wanted  to  give  an  initial  rotation  of  5°  for  the  undeformed  system,  one 
cannot  simply  set  N1R1  =  5°  because  this  would  give  an  initial  rotation  to  the  first 
node  without  propagating  this  along  the  appendage.  This  initial  rotation  will  translate 
node  2  which  will  translate  node  3  and  so  on.  Instead,  the  initial  rotation  is  done  as 
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Figure  3.2:  Simulink  Block  Diagram  For  Closed  Loop  Response 


follows: 
for  initial  roll 


NaTl  =  0 

NaT2  =  —  Alsin(NlRl) 

NaT3  =  Alcos(NlRl)  —  A  l 
NaRl  =  N1R1 
NaR2  =  0 

NaR3  =  0  (3.52) 

for  initial  pitch 

NaTl  =  Alsin{NlR2) 

NaT2  =  0 

NaT3  =  Alcos(NlR2)  —  A  l 
NaRl  =  0 
NaR2  =  NIR2 


NaR3 


0 


(3.53) 
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and  for  initial  yaw 

NaTl  =  0 

NaT2  =  0 

NaT3  =  0 

NaRl  =  0 

NaR2  =  0 

NaR3  =  N1R3  (3.54) 

which  is  similar  to  the  rotational  kinematics  presented  in  Appendix  C  where  Na  refers 
to  an  incremental  node  number  with  a  length  A l  from  the  first  node. 

Since  the  roll  and  pitch  axes  will  require  more  control  effort,  an  initial  angular 
displacement  of  5°  is  applied  to  the  pitch  axis  and  the  controller  is  trying  to  bring  the 
system  to  a  desired  attitude  of  roll  =  0°,  pitch  =  0°,  and  yaw  =  0°.  The  simulation 
generates  both  the  control  effort  and  dynamic  response  for  the  flexible  system. 

The  control  effort,  when  the  initial  design  parameters  are  used,  is  shown  in  Fig¬ 
ure  3.3.  The  peak  value  of  9.77Nm  exceeds  the  limit  of  3Nm.  This  shows  the  controller 
is  weighting  the  states  too  heavily  in  the  LQG  process  while  not  placing  enough  empha¬ 
sis  on  the  control  effort.  Since  the  feedback  gains  are  a  function  of  the  ratio  between 
Q  and  R,  one  can  either  decrease  Q  or  increase  R.  After  an  iterative  modification  to 
Q,  the  final  design  parameters  listed  in  Table  3.3  produced  the  control  effort  shown  in 
Figure  3.4. 
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Figure  3.3:  Control  Effort  with  Initial  Design  Parameters 


Table  3.3:  Final  Values  for  Design  Parameters 


Design  Parameter 

Value 

Q 

0.01  *  Inxn 

R 

hx3 

Qo 

variance  of  1 

Ro 

variance  of  0.01 

q 

[1  le2  le4  le6] 

Max  Control  Effort 

3Nm 
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Figure  3.4:  Control  Effort  with  Final  Design  Parameters 
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Stability  Robustness  Test  Based  on  the  Recovered  Compensator 


Figure  3.5:  Baseline  SR  Test  for  Initial  Design  Model 


The  stability  robustness  test,  outlined  in  Section  3.3.1,  is  evaluated  to  determine 
which  resonant  inodes  listed  in  Table  2.7  need  to  be  included  in  the  design  model.  The 
stability  robustness  (SR)  barrier  is  violated  if  only  the  first  six  modes  are  retained  in 
the  design  model  while  all  of  the  resonant  modes  are  placed  in  the  uncertainty  transfer 
matrix  (see  Figure  3.5).  Do  not  confuse  this  design  model  with  a  rigid  body  model 
since  the  resonance  frequencies  are  in  the  10~ '  to  10-5  Hz  range  and  are  not  zero.  This 
violation  of  the  SR  barrier  occurs  at  the  same  frequency  as  the  first  bending  mode. 

To  ensure  the  baseline  system  meets  SR  requirements,  additional  resonant  modes 
need  to  be  included  in  the  design  model.  Figure  3.6  is  a  plot  of  the  SR  barrier  when 
the  first  bending  and  torsional  modes  are  included.  This  design  model  has  a  stability 
margin  of  41  dB  in  the  frequency  region  of  the  second  bending  mode.  The  upper  curve 
sloping  upwards  indicates  good  tolerance  of  high-frequency  uncertainty. 

What  is  gained  by  including  the  second  bending  mode?  Stability  robustness  can 
be  improved  if  more  resonant  modes  are  included  in  the  design  model.  Performance 
of  the  closed  loop  system  depends  on  the  low-frequency  gain  and  crossover  frequency 
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Stability  Robustness  Test  Based  on  the  Recovered  Compensator 


Figure  3.6:  Baseline  SR  Test  with  First  Bending  and  Torsional  Modes 


of  the  loop  transfer  matrix  GpGc.  Larger  values  indicate  better  tracking  performance. 
Any  increase  in  gain  to  improve  tracking  performance  will  decrease  the  stability  robust¬ 
ness.  The  baseline  model  already  takes  into  consideration  the  maximum  control  effort 
produced  by  the  reaction  wheels  and  an  increase  in  bode  gain  will  exceed  this  limit. 

In  addition  to  the  system  not  being  able  to  take  advantage  of  an  increase  in  control 
gain,  the  cost  of  needlessly  adding  higher  frequency  resonant  modes  may  prove  detri¬ 
mental  to  onboard  processor  limitations.  The  LQG/LTR  process  generates  a  controller 
of  the  same  order  as  the  system  plant.  The  controller  for  a  design  model  containing  nine 
resonant  modes  is  a  3  by  3  matrix  of  transfer  functions,  each  the  ratio  of  a  17th  order 
polynomial  to  a  18th  order  polynomial.  It  would  require  programming  333  coefficients 
to  implement  this  controller.  By  adding  the  second  bending  modes  to  the  design  model, 
the  order  of  the  controller  is  increased  to  22  and  requires  an  additional  72  coefficients. 
The  inclusion  of  the  first  bending  mode  in  both  the  x  and  y  directions  and  the  first 
torsional  mode  is  enough  to  satisfy  the  stability  robustness  of  the  controller  and  the 
inclusion  of  higher  resonant  modes  is  not  required. 
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Figure  3.7:  Pole-Zero  Map  of  Each  Element  of  the  3x3  Controller  Transfer  Function 


Using  the  controller  reduction  technique  outlined  in  Section  3.3.2,  the  number  of 
required  coefficients  is  greatly  reduced.  To  illustrate  the  pole-zero  cancelation  technique, 
consider  the  pole-zero  map  of  the  controller  transfer  function  shown  in  Figure  3.7.  It  is 
apparent  several  controller  poles  lie  in  close  proximity  to  controller  zeros.  Cancelations 
may  occur  as  long  as  the  Bode  plots  of  the  reduced  controllers,  as  well  as  the  system 
dynamic  response,  match  those  of  the  unreduced  controllers.  The  Bode  plots  of  the 
unreduced  controller  are  shown  in  Figure  3.8. 

Recall  that  the  default  tolerance  value  for  the  Matlab  minreal  command  is  y/eps. 
It  is  possible  to  increase  this  tolerance  and  generate  the  same  Bode  plot  of  the  reduced 
controller  as  the  unreduced  one.  As  long  as  the  plots  for  both  the  full  and  reduced  order 
controllers  match  up,  then  the  dynamics  of  both  controllers  are  similar.  For  example, 
consider  the  controllers  if  the  default  tolerance  setting  is  used,  y/eps,  for  the  transfer 
function  going  from  input  1  to  output  1,  Gc(l,l).  From  Figure  3.9,  one  can  see  how  the 
reduced  order  controller  (green  x’s)  matches  up  with  the  full  order  controller  (red  line) 
in  both  magnitude  and  phase,  for  both  low  and  high  frequencies. 

One  may  relax  the  tolerance  in  the  minimum  realization  until  a  divergence  in  the 
Bode  plots  is  noticed.  A  better  approach  is  to  run  a  sweep  of  tolerance  values  and, 


Magnitude  (dB)  ;  Phase  (deg) 

To:Out(3)  To:Out(3)  To:Out(2)  To:Out(2)  To:Out(1)  To:Out(1) 
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Figure  3.8:  Bode  Plots  of  Each  Element  of  the  3x3  Controller  Transfer  Function 


Figure  3.9:  Comparison  of  Reduced  and  Full  Order  Controller  (default  tol) 
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Figure  3.10:  RMS  of  Bode  Magnitude  for  Gc(l,l) 


instead  of  visually  comparing  the  Bode  plots,  calculate  the  RMS  values  and  plot  those 
values  vs.  the  tolerance  step.  The  tolerance  is  stepped  from  y/eps  to  le-1  in  orders  of 
magnitude  such  that  tolerance  step  1  is  y/eps,  tolerance  step  2  is  le-7,  tolerance  step  3 
is  le-6,  and  so  on  out  to  tolerance  step  8  equalling  le-1.  The  RMS  of  the  Bode  plots  are 
calculated  for  both  magnitude  and  phase  for  all  9  controller  transfer  functions.  Since 
the  previous  plot  was  for  Gc(l,l),  the  following  two  plots  are  the  RMS  plots  for  the 
same  transfer  function. 

Looking  at  Figure  3.10  and  Figure  3.11,  a  noticeable  change  in  RMS  values  occurs 
between  tolerance  steps  3  and  4  (le-6  and  le-5).  This  means  a  noticeable  change  in 
the  reduced  order  controller  dynamics  occurs  when  additional  pole-zero  cancelations 
are  performed  beyond  a  tolerance  of  le-6.  Figure  3.12  and  Figure  3.13  show  the  Bode 
comparisons  for  both  tolerance  steps. 

Notice  how  the  two  Bode  plots  differ  once  the  RMS  is  noticeable.  This  would 
lead  one  to  believe  that  the  tolerance  limit  should  be  set  at  le-6.  However,  note  how  the 
plot  in  Figure  3.13  diverges  at  low  frequencies.  This  is  in  a  location  where  the  poles  and 
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Figure  3.11:  RMS  of  Bode  Phase  for  Gc(l,l) 


Figure  3.12:  Comparison  of  Reduced  and  Full  Order  Controller  (tol=le-6) 


102 


Figure  3.13:  Comparison  of  Reduced  and  Full  Order  Controller  (tol=le-5) 


zeros  are  in  close  proximity  to  each  other  and  the  origin  of  the  s-plane.  It  is  expected 
to  see  more  cancelations  occurring  in  this  region  than  in  the  regions  where  the  distance 
from  the  s-plane  origin  (i.e.  frequency)  is  greater. 

Looking  at  the  simulation  results  of  how  the  system  responds  to  the  reduced 
order  controller  shows  no  noticeable  difference  between  tolerance  step  3  or  4.  Instead 
of  visually  determining  this,  RMS  values  can  be  calculated  and  plotted  in  the  same 
manner  as  above  for  both  control  effort  and  dynamic  response  of  the  system.  The 
following  figures  plot  these  RMS  values. 

From  Figure  3.14  and  Figure  3.15,  a  noticeable  change  in  RMS  values  occurs 
between  tolerance  steps  6  and  7  (le-3  and  le-2).  A  design  consideration  is  to  limit 
the  pole-zero  cancelation  by  using  either  the  controller  dynamics  or  the  overall  system 
dynamics.  With  the  controller  approach,  the  reduced  controller  closely  resembles  the 
dynamics  of  the  full  order  controller.  However,  the  cost  of  this  approach  is  an  increased 
number  of  coefficients  required  to  code  up  the  controller.  If  the  system  dynamics  ap¬ 
proach  is  used,  a  smaller  number  of  coefficients  will  be  needed. 
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RMS  of  Full/Reduced  Order  Control  Effort 


Figure  3.14:  RMS  of  Full/Reduced  Order  Control  Effort 


RMS  of  Full/Reduced  Order  Dynamic  Response 


Figure  3.15:  RMS  of  Full/Reduced  Order  Dynamic  Response 
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Figure  3.16:  Stability  Robustness  Test  Using  Reduced  Controller 


To  illustrate  this,  the  design  model  is  comprised  of  the  first  nine  resonant  modes. 
For  the  full  order  system,  333  coefficients  are  required.  For  the  Bode  dynamics  approach 
to  pole-zero  cancelation,  245  coefficients  are  required.  The  system  dynamic  approach 
only  needs  108  coefficients  (67.6%  reduction  in  the  number  of  coefficients  when  compared 
to  the  full  order  system). 

The  SR  plot  using  the  reduced  controller  is  shown  in  Figure  3.16  and  generates 
the  same  stability  robustness  as  that  shown  in  Figure  3.6  for  the  full  order  controller. 
Figure  3.17  and  Figure  3.18  compare  the  full  order  controller  performance  (dashed  line) 
to  the  reduced  order  controllers  (marked  with  x’s).  From  the  singular  value  plots,  the 
control  effort,  and  dynamic  response,  it  is  apparent  that  the  reduced  order  controller 
has  the  same  dynamic  response  as  the  full  order  LQG/LTR  controller. 

Transient  system  performance  is  often  described  in  terms  of  the  unit  step  function 
response.  Since  the  input  into  the  system  is  optimally  shaped,  the  system  performance 
specifications  are  similar  with  the  following  definitions.  Dynamic  delay  time  is  the  time 
required  for  the  system  to  reach  50%  the  initial  displacement  value.  Dynamic  rise  time 
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Figure  3.17:  Control  Effort  of  Full  and  Reduced  Order  Controller 


Figure  3.18:  Dynamic  Response  of  Full  and  Reduced  Order  Controller 
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Figure  3.19:  Singular  Values  of  Recovered  Loop  Transfer  Matrix,  GcGp 


is  the  time  required  for  the  response  to  go  from  90%  to  10%  of  its  initial  value.  The  time 
which  the  system  response  settles  within  2%  of  its  initial  value  is  the  dynamic  settling 
time.  Dynamic  peak  overshoot  is  the  maximum  difference  between  the  transient  and 
steady  state  solution  and  is  a  represented  as  a  percentage  of  the  initial  displacement. 

Another  performance  specification  which  is  measured  is  the  crossover  frequency, 
uj o,  of  the  loop  transfer  matrix.  While  the  stability  margin  is  an  indicator  of  the  system 
performance,  the  crossover  frequency  determines  the  speed  of  the  system  response.  A 
higher  value  for  ujq  means  faster  response.  The  crossover  frequency  is  determined  by 
the  frequency  at  which  the  minimum  singular  value  of  the  loop  transfer  matrix,  GcGp , 
has  a  gain  of  0  dB  (see  Figure  3.19). 

System  performance  specifications  for  the  baseline  model  are  listed  in  Table  3.4. 
The  following  chapter  illustrates  how  these  performance  specifications  are  affected  when 
the  elastic  properties  of  the  appendage  and  its  length  are  modified. 
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Table  3.4:  Baseline  System  Performance 


Dynamic  Performance 

Value 

Delay  Time 

4.1  sec 

Rise  Time 

8.2  sec 

Settling  Time 

48.8  sec 

Peak  Overshoot 

21.3% 

Stability  Robustness 

41  dB 

Crossover  Frequency 

0.025  Hz 

Chapter  4 


Variation  of  Appendage  Parameters 


4.1  Satellite  Design  Concerns 

A  survey  was  presented  to  lead  satellite  design  engineers  involved  with  the  pro¬ 
grams  at  the  Space  Systems  Research  Center  (SSRC)  located  at  the  US  Air  Force 
Academy.  The  involvement  characterization  of  those  engineers  who  responded  to  the 
survey  includes  eight  satellite  control  engineers  working  directly  for  the  SSRC  and  the 
Department  of  Astronautics,  three  engineers  who  formally  worked  on  SSRC  satellite 
design  projects,  two  engineers  working  for  the  Air  Force  Research  Laboratory  in  the 
Space  Vehicles  and  Propulsion  Directorates,  two  former  SSRC  members  who  are  cur¬ 
rently  working  on  PhD  programs  with  Surrey,  and  four  first-degree  Academy  Cadets 
(undergraduate  seniors)  taking  the  small  satellite  design  course. 

The  survey  asked  each  individual  what  were  their  top  three  small  satellite  design 
concerns  from  the  following  list: 

Control  Complexity:  Number  of  resonant  modes  included  in  the  nominal  model  and 
design  model,  reduced  order  of  the  controllers,  etc. 

Material  Uncertainties:  Variations  in  elastic  modulus,  thickness,  laminate  structure, 
on-orbit  fatigue  cycles,  etc. 

The  ratio  of  the  first  bending  mode  frequency  to  the  crossover  frequency  varied  via 
length  variation  of  the  appendage. 
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Control  Design  Parameters:  Weight  on  states  and  control  effort,  loop  transfer  re¬ 
covery  gain,  diagonal  matrices,  frequency  shaped,  tuning  parameters  in  Kalman 
filter,  etc. 

Spacecraft  Parameters:  Inertia  and  mass  of  spacecraft  and  tip  mass,  appendage 
properties,  structural  damping  ratio,  appendage  configuration,  etc. 

Sensors  and  Actuators:  Differing  types  and  locations 

The  three  areas  of  concern  for  satellite  designers  surveyed  were  control  complexity 
in  the  implementation  sense,  system  response  and  sensitivity  in  the  face  of  material 
uncertainties,  and  how  the  system  performs  as  the  ratio  ^  varies  as  a  result  of  increasing 
the  length  of  the  appendage.  Small  flexible  spacecraft  performance  characteristics  used 
to  analyze  these  three  concerns  were  settling  time,  delay  time,  rise  time,  peak  overshoot, 
number  of  processor  operations  for  each  controller,  crossover  frequency  and  stability 
robustness.  A  discussion  on  each  of  the  three  primary  satellite  design  concerns  are 
presented  in  the  following  sections. 

4.2  Control  Complexity 

The  complexity  of  the  controller  does  not  refer  to  whether  a  simple  classical 
PID  controller  or  compensator  is  used  instead  of  a  more  involved  optimal  controller. 
The  calculation  of  the  controller  is  done  on  the  ground  where  computer  processing 
efforts  are  of  little  concern.  The  controllers  are  determined  and  the  coefficients  of  the 
controller  transfer  functions  are  loaded  onto  the  spacecraft’s  on-board  processor.  A 
satellite  designer’s  concern  is  if  the  on-board  processor  will  be  able  to  run  the  code  or 
not.  A  designer  would  like  to  see  as  few  lines  of  code  and  to  minimize  the  number 
of  operations  per  computing  cycle  as  possible  as  long  as  the  system  performs  within 
mission  requirements.  This  allows  valuable  processor  time  to  be  dedicated  to  other 
system  tasks  such  as  telemetry  and  data  handling,  payload  management,  and  health  and 
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status  updates.  This  section  considers  what  capabilities  a  satellite  designer  gains  from 
using  a  higher  order  controller  while  paying  increased  computing  efforts  to  implement 
such  controllers. 

The  primary  concern  is  to  measure  increased  processor  costs.  Will  a  satellite 
designer  be  able  to  use  a  current  processor  to  successfully  meet  mission  requirements 
or  will  a  more  capable  processor  be  required?  The  answer  to  this  question  lies  in 
determining  the  number  of  floating  point  operations  a  processor  can  dedicate  to  the 
attitude  control  system.  This  value  is  dependant  on  processor  memory  size  and  speed, 
other  system  functions  the  processor  needs  to  support  and  the  sequencing  of  those 
functions,  and  the  number  of  operations  per  computing  cycle  required  to  implement  the 
controller.  The  latter  processor  dependency  is  addressed  here  while  the  first  two  reside 
in  the  realm  of  system  integration  and  concept  of  operations. 

The  FLOPS  command  in  Matlab  was  used  by  satellite  designers  to  count  the 
number  of  floating  point  operations  certain  programs  and  functions  used  during  their 
execution.  However,  since  the  release  of  Matlab  version  6.0,  this  command  became 
obsolete  and  is  no  longer  practical.  In  addition,  satellite  processors  typically  run  C+ 
code  and  not  Matlab.  Time  limitations  of  the  research  efforts  also  makes  it  unfeasible 
to  benchmark  simulations  of  the  controller  on  various  commercial-off-the-shelf  (COTS) 
satellite  processors  and  operating  systems. 

The  PROFILE  command  in  Matlab  is  useful  in  determining  the  execution  time  of 
a  program  or  function.  It  is  mainly  used  to  debug  and  optimize  run  times  of  M-files  by 
providing  information  such  as  execution  time,  number  of  calls,  parent  functions,  child 
functions,  code  line  hit  count,  and  code  line  execution  time.  While  useful  when  running 
Matlab  commands,  the  PROFILE  command  does  not  track  Simulink  simulation  runs. 

The  TIC  and  TOC  commands  can  be  used  while  executing  a  Simulink  model  from 
the  Matlab  command  window.  Tic/toc  are  stopwatch  commands  which  provide  elapsed 
time  measurements  from  when  tic  starts  the  watch  and  when  toe  stops  it.  The  command 
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entered  in  the  Matlab  command  window  is  tic,[t,x,y]=sim(’modelname’,runtime);toc. 
This  command  starts  the  stopwatch,  executes  the  Simulink  model  ’modelname’  and 
runs  from  0  to  runtime  (in  seconds),  then  stops  the  stopwatch  and  prints  the  elapsed 
time. 

Test  runs  of  the  Simulink  model  were  for  60  second  simulations.  Since  integrators 
like  ode45  have  variable  time  steps,  it  is  necessary  to  change  the  time  step  size  in  the 
simulation  configuration  parameters  solver  options  from  variable-step  sizes  to  fixed  time 
steps  so  that  each  run  was  conducted  over  fixed  time  intervals.  The  Runge-Kutta  solver 
ode4  was  used  with  fixed  time  step  size  of  1  ms.  In  addition,  various  programs  open  on 
the  desktop  may  use  CPU  resources  while  the  simulations  were  being  run.  Therefore, 
prior  to  data  collection,  the  computer  needs  to  be  restarted  with  only  Matlab  open.  This 
will  prevent  occasional  pings  to  the  system  from  the  network  or  resource  allocations  used 
from  Outlook  during  received  emails. 

Even  after  restarting  the  computer,  slight  variations  in  run  time  may  be  experi¬ 
enced  during  several  runs  of  the  same  simulation.  Larger  variations  were  noticed  when 
the  CPU  usage  history  demonstrated  increased  numbers  of  non-linear  affects  instead  of 
smooth  histories.  An  interesting  observation  is  larger  occurrences  of  non-linear  CPU 
usage  was  noted  when  simulations  involved  larger  numbers  of  controller  coefficients  (see 
Figure  4.1  and  Figure  4.2).  This  is  meaningful  in  that  if  the  on-board  processor  must 
cache  data,  a  non-linear  usage  of  the  processor  will  be  noted. 


Figure  4.1:  Smooth  CPU  Usage  History 


To  determine  a  more  accurate  estimate  of  run  times,  10  runs  were  done  for  each 
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simulation  with  the  average  time  being  recorded.  In  addition,  each  run  wasn’t  started 
until  the  hard  drive  was  no  longer  being  accessed  and  the  CPU  Usage  History  was  0% 
for  5  seconds. 

Elapsed  execution  times  were  recorded  for  differing  design  models  based  on  in¬ 
cluded  flexible  modes,  m.  From  previous  work  presented  in  Chapter  3,  the  truth  model 
includes  the  first  10  flexible  modes,  m=10,  along  with  the  six  low  frequency  resonant 
modes  while  the  baseline  design  model  included  the  first  flexible  modes  about  each  axis, 
m=3.  Not  only  were  simulations  run  for  the  full  order  controllers  for  each  design  model, 
but  times  were  also  calculated  for  reduced  order  controllers. 

The  number  of  operations  performed  per  computing  cycle  is  equivalent  to  the 
number  of  addition,  subtraction,  multiplication,  and  division  operations  conducted.  If 
a  transfer  function  of  the  controller  is  of  the  general  form 

clqs  T  ais  T  ■  ■  ■  T  o,m— is  +  am  .  . 

bosn  T  b\sn  1  +  •  •  •  +  bn-is1  +  bn 

where  a  and  b  are  coefficients  needing  to  be  programmed,  n  is  the  order  of  the  controller, 
and  m  is  the  order  of  zeros.  Recall  there  are  nine  controller  transfer  functions  for 
a  3-input  3-output  system.  Considering  only  the  addition  and  multiplication  of  the 
coefficients,  the  total  number  of  operations  performed  is  equal  to 

9  9 

operations  =  E  2(ni)  +  ^2(mi)  +  18  (4.2) 

i  i 

Since  m  =  n  —  1 ,  and  adding  to  the  result  nine  for  the  division  within  the  transfer 
functions,  six  for  the  computation  of  each  of  the  three  input  signal  as  a  result  of  three 


Figure  4.2:  Non-Linear  CPU  Usage  History 
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measurement  contributions  for  each,  and  adding  an  additional  three  for  feeding  back 
the  three  inputs,  Eq.4.2  becomes 

9  9 

operations  =  E  2 (nj)  +  2 {rii  —  1)  +  18  +  9  +  6  +  3 

i  i 

9 

=  ^4(ni)  +  18  (4.3) 


As  was  shown  in  Chapter  3,  the  number  of  coefficients  needing  to  be  programmed 
into  the  on-board  processor  is 

9 

coefficients  =  ^  2  (m)  +  9  (4.4) 

i 


and  the  number  of  operations  performed  per  computing  cycle  is  twice  the  number  of 
coefficients  of  the  controller  transfer  function  matrix. 

The  elapsed  run  time  for  each  case  is  plotted  vs  the  number  of  operations  for 
each  design  model  and  is  shown  in  Figure  4.3.  The  legend  identifies  which  plot  line 
corresponds  to  which  design  model  where  m  is  the  value  representing  the  number  of 
flexible  modes  included  in  the  design  model.  The  plot  is  linear  except  in  the  region 
where  the  number  of  operations  exceeds  600.  This  is  a  result  of  the  non-linear  usage  of 
the  CPU  capability  and  caching  of  virtual  memory  to  the  hard  drive. 

The  data  collected  for  the  run  times  was  collected  on  a  Pentium(R)  4  CPU,  1.69 
GHz,  with  256  MB  of  RAM  and  using  the  Microsoft  Windows  XP  version  2002  opera¬ 
tion  system.  However,  satellite  processors  are  not  using  the  same  processor  and  will  run 
code  in  forms  other  than  M-files.  To  account  for  this,  the  data  is  non-dimensionalized 
by  dividing  through  all  of  the  run  times  by  the  time  it  took  to  execute  the  baseline  de¬ 
sign  model  with  the  smallest  order  controller  (pole-zero  cancellations  determined  with 
tolerance  equal  to  0.1  during  the  minimum  realization  calculations).  This  allows  a  com¬ 
parison  of  implementation  difficulty  of  one  design  model/reduced  controller  over  another 
(see  Figure  4.4).  This  information  is  useful  for  satellite  designers  to  determine  on-board 
processor  size  and  capabilities  while  considering  more  complex  control  structures  by 
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Figure  4.3:  Elapsed  Time  vs  Number  of  Operations 
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Figure  4.4:  Implementation  Costs  of  Various  Order  Design  Models 


allowing  them  to  quickly  look  up  how  many  operations  are  required  to  implement  the 
various  orders  of  controllers. 

Figure  4.5  takes  a  look  at  how  the  number  of  operations  per  computing  cycle 
decreases  as  the  tolerance  step  changes.  Using  this  plot,  satellite  designers  can  consider 
how  many  operations  their  on  board  processor  can  dedicate  to  the  attitude  control 
architecture  and  determine  which  design  model  to  implement  and  how  far  they  need  to 
reduce  the  controller  order  before  needing  to  acquire  a  higher  performing  processor.  For 
example,  if  a  satellite  designer  determines  that  a  given  processor  is  capable  of  dedicating 
500  operations  per  cycle  to  the  control  system,  then  a  design  model  which  includes  the 
first  three  flexible  modes  may  be  run  regardless  of  the  reduced  order  of  the  controller. 
However,  to  implement  any  of  the  other  design  models  which  include  more  flexible 
modes,  the  controllers  need  to  be  reduced  in  order  while  using  tolerance  equal  to  0.0001 
during  the  minimum  realization  process.  Now,  the  satellite  designer  must  determine 
whether  they  need  to  go  to  a  better  performing  processor  or  not. 

The  main  considerations  for  a  satellite  designer,  when  determining  processor  size, 
is  will  the  added  cost  and  power  consumption  generate  increased  stability  robustness 
and  performance.  As  shown  in  Chapter  3,  the  baseline  model  (m=3)  had  a  stability 
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Figure  4.5:  Pole-Zero  Cancelation  Effects  on  Number  of  Operations 
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Figure  4.6:  Dynamic  Performance  and  Implementation  Difficulty 


robustness  of  41  dB.  Adding  more  flexible  modes  in  the  design  model  increases  this 
stability  robustness.  However,  since  the  reaction  wheels  have  a  limit  placed  on  them  with 
how  much  torque  they  can  generate,  the  increased  stability  robustness  gains  the  satellite 
designer  nothing  since  the  bode  gain  cannot  be  increased  to  improve  performance  (the 
control  effort  limit  of  3Nm  and  how  the  weighting  matrices  in  the  LQG/LTR  design 
impact  this  effort  and  system  performance  is  shown  in  Chapter  3). 

Since  a  satellite  designer  doesn’t  benefit  from  an  increase  in  stability  robustness, 
the  consideration  of  dynamic  performance  becomes  important.  Now,  the  full  order 
controller  generated  for  the  truth  model  is  the  benchmark  when  determining  system 
performance.  RMS  values  are  calculated  for  each  design  model  and  each  reduced  con¬ 
troller  to  compare  their  dynamic  response  to  that  of  the  full  order  truth  model.  A 
3-dimensional  plot  is  created  to  show  the  relationship  between  these  RMS  values,  im¬ 
plementation  difficulty,  and  number  of  operations  (see  Figure  4.6). 

A  satellite  designer  can  use  Figure  4.6  to  determine  which  control  complexity  to 
use  on  the  satellite.  A  favorable  location  occurs  in  the  lower/middle  portion  of  the 
graph  where  RMS  values  of  the  reduced  system  response  compared  to  the  full  order 
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Table  4.1:  Performance  Comparison  Between  Baseline  and  Truth  Model 


Performance 

Baseline  Model 

Truth  Model 

Delay  Time  (sec) 

4.1 

4.1 

Rise  Time  (sec) 

8.2 

8.3 

Settling  Time  (sec) 

48.8 

48.8 

Peak  Overshoot  (%) 

21.3 

21.0 

Crossover  Frequency  (Hz) 

0.025 

0.025 

Stability  Robustness  (dB) 

41 

125 

truth  model  response  are  the  smallest  with  requiring  fewer  operations  per  cycle  and 
easier  implementation.  From  this  plot,  and  knowing  that  the  baseline  model  (m=3) 
meets  stability  robustness  requirements,  the  baseline  model  presented  in  Chapter  3  is 
an  ideal  choice  for  a  satellite  designer.  The  system  performance  of  the  baseline  model 
is  nearly  identical  to  the  performance  of  the  full  order  truth  model  except  for  stability 
robustness  (see  Table  4.1).  The  same  performance  is  achieved  with  the  baseline  model 
while  using  37.2%  the  number  of  operations,  being  20  times  easier  to  implement,  and 
only  has  an  dynamic  RMS  error  of  0.01°.  Therefore,  the  baseline  model  is  the  better 
choice  unless  mission  requirements  desire  pointing  accuracies  of  0.01°  or  better. 

4.3  Material  Uncertainty 

System  performance  was  nearly  identical  between  the  baseline  model  using  re¬ 
duced  order  controllers  and  the  truth  model  using  full  order  controllers.  Satellite  de¬ 
signers  can  use  this  information  in  the  sizing  of  on-board  processors.  However,  another 
concern  the  survey  indicated  was  would  mission  requirements  be  greatly  impaired  if 
material  uncertainties  in  the  flexible  structure  were  taken  into  consideration. 

Not  every  elastic  memory  composite  component  is  manufactured  with  identical 
properties.  A  certain  level  of  uncertainty,  or  error  bar,  exists  between  one  piece  and 
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another.  CTD  mentions  this  is  a  result  of  variances  in  material  thickness,  laminate  prop¬ 
erties,  mesh  formulation,  variances  in  parasitic  components  such  as  heating  coils,  resin, 
connections,  etc.  The  composite  industry  typically  does  not  measure  every  possible 
variance  in  the  material.  Instead,  they  lump  the  material  uncertainties  into  the  elastic 
modulus  and  accept  an  uncertainty  of  ±10%  of  the  baseline  elastic  modulus  (3-sigma). 

One  consideration  which  is  not  included  in  the  error  bar  is  on-orbit  fatigue  cycling 
of  the  composite  material.  Currently,  no  data  exists  (theoretical  or  experimental)  to 
determine  a  level  of  uncertainty  in  the  elastic  modulus  resulting  from  prolonged  exposure 
to  the  space  environment.  On-orbit  fatigue  of  the  material  may  result  in  an  increase  of 
stiffness  (non-linear  stress  to  strain  relationships,  resin  hardening,  etc)  or  an  increase  in 
flexibility  (severed  laminate  fibers,  variance  in  thermal  transition  barriers,  etc)  in  the 
material. 

To  account  for  the  unmodeled  fatigue  cycles,  numerical  analysis  of  the  baseline 
model  was  run  while  varying  the  baseline  elastic  modulus,  E0,  from  80%  to  120%  at  1% 
intervals.  This  allowed  the  modified  elastic  modulus,  E) ,  to  take  on  values  of  0.80Eo, 
0.81Eo,  1.19E0,  and  1.20Eo.  While  the  ratio,  was  varied,  the  following  para¬ 

meters  were  held  constant:  mass,  inertia,  structural  shape,  area,  poisson  ratio,  shear 
modulus,  density,  weighting  matrices  Q  and  R,  recovery  gain  r,  and  the  tolerance  used 
in  the  minimum  realization  of  the  controllers. 

An  analytical  relationship  of  how  the  resonant  frequencies  change  as  a  function 
of  elastic  modulus  is  shown  in  the  following  derivation: 

[ea  ,  s 

(4-5) 

where  u  is  the  baseline  resonant  frequency,  c  is  a  characteristic  value  dependent  on 
boundary  conditions  and  resonant  number,  E  is  the  elastic  modulus,  A  is  the  cross- 
sectional  area,  m  is  the  appendage  mass,  and  L  is  its  length. 

When  the  elastic  modulus  is  varied  by  a  value  of  SE,  the  modified  resonant 
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frequency,  co*,  is  found  by 


co  =  c\ 
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(4.6) 


(4.7) 


Analytical  relationships  between  elastic  modulus  variation  and  system  perfor¬ 
mance  characteristics  do  not  exist.  However,  it  is  possible  to  estimate  these  relation¬ 
ships  by  taking  data  from  the  numerical  analysis  on  the  baseline  model  as  perviously 
mentioned.  Each  elastic  modulus  ratio  is  entered  into  Patran  and  generates  a  new 
eigen  analysis.  This  information  is  placed  into  an  input  file  for  the  Matlab  code  and 
dynamic  response  simulations  are  executed.  Values  for  peak  control  effort,  stability  ro¬ 
bustness,  crossover  frequency,  delay  time,  rise  time,  settling  time,  and  percent  overshoot 
were  recorded.  The  data  was  plotted  versus  the  elastic  modulus  ratios  to  determine  a 
relationship  between  material  uncertainties  and  system  performance. 

The  Basic  Fitting  Option  was  used  to  evaluate  the  plots  and  generate  a  ’’best  fit” 
polynomial,  in  the  least  squares  sense,  for  a  given  set  of  data.  Residuals  were  calculated 
as  a  measure  of  how  well  the  predicted  data  matches  the  observed  data.  If  the  residuals 
show  strongly  patterned  behavior,  then  it  should  be  possible  to  do  better  than  a  simple 
polynomial  fit  (exponential  fit).  Error  bounds  are  also  calculated  for  each  data  set 
to  determine  if  the  data  is  reasonably  modeled  by  the  fit.  The  error  bound  uses  an 
interval  of  ±2 6  which  corresponds  to  a  95%  confidence  interval.  The  Matlab  commands 
POLYFIT  and  POLYVAL  were  used  to  perform  curve  fitting. 

As  the  flexible  structure  increases  in  stiffness,  the  crossover  frequency  of  the  loop 
transfer  matrix  remained  constant,  peak  control  effort,  stability  robustness,  and  peak 
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Table  4.2:  Performance  Characteristics  as  a  Function  of  Material  Uncertainty 


Performance  Characteristic 

Material  Uncertainty  Ratio  (x) 

Peak  Control  Effort  (Nm) 

=  1.3551  +  1.6588x 

Stability  Robustness  (dB) 

=  31.136  +  8.1023x 

Crossover  Frequency  (Hz) 

=  0.025 

Delay  Time  (sec) 

=  4.5491  -  0.46261x 

Rise  Time  (sec) 

=  9.2476  -  0.9941x 

Settling  Time  (sec) 

=  49.729  -  1.0029x 

Peak  Overshoot  (%) 

=  17.053  +  4.4043x 

overshoot  increases,  and  delay,  rise,  and  settling  times  decrease.  The  linear  equations 
for  these  relationships  are  listed  in  Table  4.2  while  the  resulting  plots  are  shown  in 
Figures  4.7  through  4.20. 
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Figure  4.11:  Linear  Fitting  of  Crossover  Frequency  Data 
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Figure  4.12:  Error  Bounds  for  Crossover  Frequency  Data 
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An  interesting  observation  is  drawn  from  the  residual  plots  in  that  each  data 
entry  resulting  from  the  Simulink  run  of  ratio  1.17  resulted  in  an  observed  value  which 
falls  outside  the  error  bound.  It  is  expected  for  some  values  to  fall  outside  the  error 
bounds  since  the  confidence  interval  is  95%.  However,  for  all  simulation  response  points 
for  this  one  case  may  indicate  a  faulty  assumption. 

A  likely  assumption  which  contributes  to  this  observation  is  holding  the  tolerance 
value  constant  during  the  minimum  realization  of  the  controllers.  The  comparison  plot 
of  the  reduced  and  full  order  controller  bode  dynamics  in  Section  3.4  shows  a  divergence 
at  low  frequencies  (see  Figure  3.13).  This  difference  is  in  the  low  frequency  region  of  the 
system  where  the  system  poles  and  zeros  are  in  close  proximity.  When  Patran  generates 
the  modal  information,  it  is  unable  to  provide  zero  frequency  information.  Numerical 
artifacts  in  the  analysis  produce  frequencies  of  1CP11  Hz.  It  is  possible  these  numerical 
artifacts  may  become  important  in  the  system  dynamic  response. 

RMS  values  were  calculated  of  both  control  effort  and  dynamic  response  for  the 
case  where  the  elastic  modulus  ratio  equaled  1.17  (see  Figure  4.21  and  Figure  4.22).  A 
large  gradient  is  apparent  in  the  tolerance  step  region  where  the  baseline  model  plots 
indicated  a  tolerance  value  of  0.001,  tolerance  step  of  6,  during  minimum  realization  of 
the  controllers  while  the  following  plots  indicate  the  tolerance  value  should  equal  10-6, 
tolerance  step  of  3. 

For  more  accurate  predictions,  the  use  of  a  full  order  controller  is  required  so  that 
low  frequency  dynamics  are  not  lost.  However,  it  turns  out  the  system  is  not  sensitive 
to  the  differing  controller  order.  The  system  performance  characteristics,  for  the  case 
where  =  1.17,  shown  in  Table  4.3  for  a  full  order  controller  and  those  of  a  controller 
realized  when  tolerance  equals  0.001  indicate  minimal  impact  of  increased  pole  zero 
cancelations  near  the  origin.  Five  of  the  seven  values  were  identical  while  the  difference 
in  settling  time  was  0.6%  and  peak  overshoot  was  0.4%. 
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Figure  4.21:  Full/Reduced  Order  Control  Effort  for  Material  Uncertainty  Ratio  1.17 


Figure  4.22:  Full/Reduced  Order  Dynamic  Response  for  Material  Uncertainty  Ratio 
1.17 
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Table  4.3:  Performance  Comparison  Between  Full  and  Reduced  Controller 


Performance 

Full  Order  Controller 

Reduced  Order  Controller 

Delay  Time  (sec) 

3.6 

3.6 

Rise  Time  (sec) 

7.1 

7.1 

Settling  Time  (sec) 

47.8 

48.1 

Peak  Overshoot  (%) 

27.4 

27.3 

Peak  Control  Effort  (Nm) 

3.65 

3.65 

Crossover  Frequency  (Hz) 

0.025 

0.025 

Stability  Robustness  (dB) 

40.3 

40.3 

4.4  Appendage  Length 

As  the  length  of  a  gravity  gradient  boom  is  increased,  the  mass  located  at  its 
deployed  end  can  be  reduced  while  maintaining  the  required  moment  of  inertia  in  the 
nadir /zenith  axis.  This  is  a  concern  for  small  satellite  designers  because  of  the  mass 
and  volume  limitations  these  satellites  possess.  While  uncertainties  in  the  appendage 
material  produced  linear  approximations  for  system  performance,  this  may  not  remain 
the  case  as  its  length  is  increased.  The  appendage  may  reach  a  certain  length  at  which 
the  design  model  will  need  to  be  modified  to  include  higher  frequency  mode  shapes. 

To  vary  the  appendage  length  for  data  collection  purposes,  a  small  step  size  (on 
the  order  of  1%  of  the  original  length)  is  used  and  system  performance  parameters 
are  calculated.  The  length  is  increased  by  the  same  step  size  and  the  parameters  are 
calculated  again.  If  the  data  appears  linear,  then  the  step  size  is  doubled.  Increasing 
the  step  size  in  this  matter  is  repeated  until  the  data  no  longer  appears  linear.  Then, 
the  initial  step  increment  is  applied  to  the  appendage  length  prior  to  the  last  step 
size  increase  and  the  process  is  continued.  The  gradient  decent  method  provides  a 
systematic  approach  to  increasing  the  appendage  length  while  spending  less  time  in 
regions  of  minimal  change  and  focusing  most  of  the  data  collection  efforts  where  the 
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data  varies  at  a  greater  rate. 

A  primary  concern  for  this  section  is  how  does  the  flexible  system  perform  as  the 
resonant  frequencies  approach  the  crossover  frequency  of  the  loop  transfer  matrix.  A 
relationship  between  the  original  resonant  frequency,  uq,  and  the  frequency  once  the 
length  is  modified,  102,  is  shown  in  the  following  equations: 
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where  the  frequency  equation  was  shown  in  Eq.  4.5. 

The  baseline  system  has  the  first  bending  mode  at  a  frequency  of  3.1518  Hz  (see 
Table  2.7)  and  a  crossover  frequency  at  0.025  Hz  (see  Table  3.4).  This  yields  a  baseline 
ratio  of  ^  =  126.07.  Inserting  this  value  into  Eq.  4.8  creates  a  method  of  relating  an 
increase  in  length  to  a  desired  frequency  ratio. 
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An  initial  run  of  the  simulations  considered  increasing  the  length  of  the  appendage 
while  holding  all  other  variables  constant  across  the  simulations.  Data  collected  is  shown 
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Table  4.4:  Simulation  Results  When  Only  Appendage  Length  is  Varied 


Length  (m) 

Peak  CE  (Nrn) 

SR  (dB) 

uj\  (Hz) 

u;0  (Hz) 

Frequency  Ratio 

4.0 

3.06 

40.7 

3.1518 

0.0241 

130.78 

4.111 

3.06 

39.6 

3.0979 

0.0250 

124.10 

4.32 

2.92 

35.4 

3.0009 

0.0242 

124.15 

4.5 

2.80 

34.8 

2.9211 

0.0234 

124.59 

4.8 

2.71 

32.8 

2.7952 

0.0228 

122.84 

5.1 

2.63 

29.6 

2.6764 

0.0215 

124.74 

5.7 

2.45 

26.0 

2.4557 

0.0211 

116.40 

6.9 

2.12 

15.3 

2.0638 

0.0181 

113.90 

7.8 

2.03 

11.8 

1.8046 

0.0176 

102.30 

8.0 

1.87 

10.7 

1.7509 

0.0164 

106.47 

8.2 

1.77 

10.9 

1.6986 

0.0164 

103.25 

in  Table  4.4.  The  heading  for  the  table  lists  appendage  length,  peak  control  effort, 
stability  robustness,  frequency  of  the  first  bending  mode,  crossover  frequency,  and  the 
ratio  of  first  bending  mode  to  crossover. 

The  stability  robustness  value  steadily  decreases  as  the  appendage  length  in¬ 
creases.  There  occurs  a  length  at  which  the  baseline  design  model  will  no  longer  remain 
robust  in  the  face  of  unmodeled/uncertain  high  frequency  structural  modes  and  noise. 
To  estimate  at  what  point  this  occurs,  the  stability  robustness  values  were  plotted  versus 
length  (see  Figure  4.23). 

If  the  decreasing  values  of  the  norm  of  the  residuals  shown  in  Table  4.5  are  used 
to  indicate  which  order  of  polynomial  to  fit  the  observed  data,  a  fifth  order  polynomial 
can  be  determined  to  estimate  at  what  length  the  SR  test  is  violated  and  should  not  be 
used  to  accurately  estimate  values  beyond  the  immediate  range  of  the  data  observed. 
The  fifth  order  polynomial  produced  is 

SR  =  —0.18411  L5  +  5.8804L4  -  73.717L3  +  453.7L2  -  1379.6L  +  1701.3  (4.10) 


Using  Eq.  4.10,  one  can  estimate  at  what  appendage  length  will  the  stability 
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Figure  4.23:  Curve  Fitting  of  Stability  Robustness  Data  for  Length  Variation 


robustness  value  be  negative,  which  indicates  the  need  to  include  the  second  bending 
modes  in  the  design  model.  This  resulted  in  an  appendage  length  of  9.2  m. 

Since  the  weighting  matrices  used  in  the  LQG/LTR  process,  Q  and  R,  were  held 
constant,  both  peak  control  effort  and  crossover  frequency  decreased.  The  intent  is  to 
maximize  the  control  input  of  the  reaction  wheels  (3.0  Nrn).  Therefore,  the  simulations 
were  repeated  while  increasing  the  q/r  ratio  so  that  the  peak  control  effort  of  3.0  Nm  is 
maintained.  The  resulting  data  is  presented  in  Table  4.6.  The  inclusion  of  the  frequency 
of  the  second  bending  modes,  u>2,  is  included  to  illustrate  how  its  ratio  changes. 

By  allowing  the  simulation  to  reach  the  control  limit  of  the  reaction  wheels,  the 
^  ratio  decreased  at  a  faster  rate.  Also,  the  stability  robustness  value  becomes  an 
issue  at  a  shorter  appendage  length  than  the  estimated  9.2m.  When  the  SR  values  of 
3.9  dB  and  3.8  dB  (for  lengths  of  8.0m  and  8.2m)  are  taken  in  consideration  with  the 
information  provided  in  Section  4.3,  a  different  conclusion  is  made.  The  linear  equation 
relating  stability  robustness  to  material  uncertainty,  the  associated  error  bounds,  and 
the  norm  of  the  residuals  indicate  the  SR  value  of  3.9  dB  for  an  appendage  length  of 
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Table  4.5:  Norm  of  Residuals  for  Appendage  Length  Variation 


Curve  Fitting  Order 

Norm  of  Residuals 

Linear 

4.7192 

Quadratic 

2.6403 

Cubic 

2.5253 

4th  order 

1.9473 

5th  order 

1.5262 

6th  order 

1.5258 

7th  order 

1.5134 

8.0m  may  not  provide  a  safe  enough  margin  for  the  baseline  design  model. 

To  further  study  when  inclusion  of  the  second  bending  modes  in  the  design  model 
is  desirable,  consider  the  plot  of  q/r  ratio  versus  appendage  length  (see  Figure  4.24).  The 
plot  is  linear  until  an  appendage  length  of  8.0  m  is  reached.  This  is  another  indication 
that  for  appendage  lengths  equal  to  or  greater  than  8.0m,  the  system  will  not  remain 
robust  in  the  face  of  unmodeled/uncertain  high  frequency  dynamics  and  noise. 

Consider  what  happens  to  the  stability  robustness  plot  for  a  design  model  using 
an  appendage  length  of  8.2m.  From  the  plots  shown  in  Figure  4.25  and  Figure  4.26,  a 
satellite  designer  sees  that  a  better  stability  robustness  value  is  achieved  when  the  second 
bending  modes  are  included  in  the  design  model.  Table  4.7  contains  key  simulation  data 
for  both  design  models  with  an  appendage  length  of  8.2m. 

When  the  second  bending  modes  are  included  in  the  design  model,  the  SR  value 
changes  from  3.8dB  to  34.7dB,  while  the  q/r  ratio  and  crossover  frequency  remain 
relatively  the  same.  Once  the  appendage  length  reaches  8.0m,  the  satellite  designer 
must  consider  stability  robustness  of  the  system  with  information  provided  on  material 
uncertainties  and  the  increased  number  of  operations  to  execute  such  a  controller  using 
on-board  processors. 
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Table  4.6:  Simulation  Results  When  Control  Effort  Limit  is  Maintained 


Length 

q/r 

Peak  CE 

SR 

w0 

Wl 

UJl_ 

WO 

w2 

W2 

WO 

4.0 

0.010 

3.06 

40.7 

0.0241 

3.1518 

130.78 

12.5924 

522.50 

4.111 

0.010 

3.06 

39.6 

0.0250 

3.0979 

124.10 

12.1313 

485.25 

4.32 

0.011 

2.99 

35.2 

0.0247 

3.0009 

121.49 

11.3213 

458.35 

4.5 

0.013 

2.99 

34.2 

0.0251 

2.9211 

116.59 

10.6764 

425.36 

4.8 

0.015 

2.99 

31.9 

0.0252 

2.7952 

110.93 

9.7202 

385.72 

5.1 

0.017 

3.00 

28.5 

0.0244 

2.6764 

109.62 

8.8757 

363.76 

5.7 

0.021 

2.99 

24.4 

0.0255 

2.4557 

96.43 

7.4989 

294.07 

6.9 

0.033 

3.00 

12.6 

0.0241 

2.0638 

85.77 

5.6494 

234.42 

7.8 

0.041 

3.00 

8.6 

0.0251 

1.8046 

71.98 

4.7640 

189.80 

8.0 

0.059 

2.99 

3.9 

0.0249 

1.7509 

70.42 

4.6073 

185.03 

8.2 

0.084 

2.99 

3.8 

0.0271 

1.6986 

62.79 

4.4592 

164.55 

Figure  4.24:  LQG/LTR  Design  Ratio,  q/r,  as  Appendage  Length  is  Varied 
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Stability  Robustness  Test  Based  on  the  Recovered  Compensator 


Figure  4.25:  SR  Plot  for  8.2m  Appendage  with  1st  Bending  Modes  in  Design  Model 


Stability  Robustness  Test  Based  on  the  Recovered  Compensator 


Figure  4.26:  SR  Plot  for  8.2m  Appendage  with  2nd  Bending  Modes  in  Design  Model 


Chapter  5 


Research  Summary  and  Recommendations 


5.1  Research  Summary 

The  trend  in  utilizing  small  satellites  to  accomplish  space  missions  has  been 
steadily  increasing  over  the  last  two  decades.  Businesses,  governmental  organizations, 
and  academic  institutions  find  the  reduced  development  costs  and  time  lines,  when 
compared  to  the  larger  conventional  satellites,  an  appealing  benefit  when  establishing 
a  small  satellite  program. 

As  the  number  of  small  satellite  missions  increase  in  the  coming  years,  so  too  will 
the  unique  ways  in  which  designers  prepare  for  these  missions.  Working  within  mass, 
volume,  and  power  constraints,  satellite  designers  will  ’’push  the  envelope”  on  what 
is  possible  to  accomplish.  These  efforts  will  generate  creative  ways  of  designing  satel¬ 
lites  to  successfully  meet  mission  requirements.  Non-traditional  methods  of  deploying 
structures  on  small  satellites  is  one  such  emerging  area  of  study. 

Current  research  efforts  conducted  within  the  materials  industry  are  looking  into 
constructing  deployable  structures  from  elastic  memory  composites  (EMC).  Strain  en¬ 
ergy  is  stored  within  appendages  made  from  shape  memory  composites  and  can  be 
released  upon  command  by  heating  the  material  beyond  its  glass  transition  tempera¬ 
ture.  Shape  memory  mechanisms  can  eliminate  the  need  for  traditional  highly  complex 
mechanical  deployment  devices,  massive  launch  canisters,  and  independent  deployment 
control  systems. 
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While  EMC  appendages  require  less  mass  and  volume  compared  to  similar  per¬ 
forming  traditional  booms  made  from  beryllium  copper,  the  increased  flexible  nature  is 
a  concern  for  the  satellite’s  attitude  control  system.  Historically,  research  in  the  area  of 
flexible  spacecraft  control  has  focused  on  assuming  cantilevered  boundary  conditions  at 
the  connection  point  between  the  spacecraft  bus  and  the  flexible  appendage.  This  as¬ 
sumption  simplifies  several  terms  in  the  Lagrangian  dynamics  because  the  center  of  mass 
of  the  system  experiences  small  variations.  However,  the  cantilever  assumption  becomes 
less  valid  as  the  mass  of  the  satellite  bus  is  reduced  and  the  appendage  configuration 
does  not  take  on  an  asymmetric  shape. 

As  the  system’s  center  of  mass  moves  further  away  from  the  center  of  mass  of  the 
controlling  body,  larger  displacements  of  the  system’s  center  of  mass  are  realized  in  the 
body  reference  frame  for  the  deformed  spacecraft.  This  means  the  time  rate  of  change 
of  the  system’s  inertia  matrix  is  not  zero  and  terms  once  neglected  in  the  Lagrangian 
now  need  to  be  taken  into  consideration.  Assuming  a  free-free  boundary  condition  for 
determining  the  mode  shapes  of  the  small  flexible  spacecraft  is  the  correct  approach  in 
creating  a  more  accurate  dynamic  model  of  the  system. 

LQG/LTR  optimal  control  techniques  were  applied  to  a  dynamic  model  generated 
from  generalized  mass  and  stiffness  matrices  created  from  a  finite  element  model  of  a 
small  satellite  using  a  gravity  gradient  boom  comprised  of  EMC  materials.  The  satel¬ 
lite  mirrors  the  FalconSat-3  spacecraft  designed  at  the  US  Air  Force  Academy’s  Space 
Systems  Research  Center.  The  design  considerations  for  the  LGQ/LTR  controller  is  to 
return  the  spacecraft  to  nominal  pointing  requirements  following  an  initial  displacement 
that  could  result  from  momentum  dumping  of  the  on-board  reaction  wheels. 

This  research  provides  insight  in  the  area  of  small  flexible  spacecraft  which  a  satel¬ 
lite  designer  can  use  to  determine  computer  processor  sizing,  how  material  uncertainties 
and  fatigue  cycling  may  impact  system  performance  parameters,  and  what  happens  to 
the  robustness  of  the  control  system  as  the  length  of  the  appendage  is  varied. 
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5.2  Future  Research  Recommendations 

Several  key  areas  of  future  research  are  mentioned  in  this  section.  A  significant 
finding  from  this  research  study  is  if  one  were  to  create  the  nominal  and  design  models 
using  the  method  outlined  in  Section  2.3,  and  then  calculate  the  additive  uncertainty 
matrix  using 

ACr  =  Gnom  Gdes  (5.1) 

then  it  becomes  necessary  to  include  the  first  bending  and  torsional  modes  in  the  design 
model  even  though  these  modes  are  two  orders  greater  than  the  magnitude  of  the 
crossover  frequency.  According  to  the  literature,  control  system  designers  model  a 
system  as  a  rigid  body  as  long  as  the  first  resonant  mode  is  at  least  one  order  of 
magnitude  greater  than  the  crossover  frequency.  Yet,  for  this  research,  the  low  frequency 
design  model  should  not  be  confused  with  a  rigid  body  model  since  the  FEM  analysis 
generates  low  frequency  information  instead  of  zero  mode  information. 

The  rule  of  thumb  of  one  order  of  magnitude  is  using  the  same  assumption  of 
boundary  conditions  as  has  been  used  repeatedly  for  the  last  40  years  of  research.  This 
assumption  is  that  the  system  is  modeled  as  a  cantilever  appendage  where  the  system 
center  of  mass  is  contained  within  the  controlling  body  and  only  experiences  a  small 
variance  in  its  location.  This  research  is  conducted  on  a  free-free  system  where  the 
center  of  mass  is  not  restricted  within  the  controlling  body  and  will  experience  larger 
variances  in  its  location.  The  one  order  of  magnitude  statement  applies  to  and  has 
been  validated  for  cantilever  systems.  Further  research  effort  is  required  to  characterize 
the  system  modeled  in  this  document  and  determine  what  factors  affect  the  stability 
robustness  of  a  free-free  system. 

One  area  of  research  currently  being  conducted  at  the  Air  Force  Research  Labs, 
Propulsion  Directorate,  is  the  application  of  micro  pulsed  plasma  thrusters  (/xPPTs) 
in  augmenting  the  attitude  control  system  of  small  flexible  spacecraft.  Future  efforts 
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could  evaluate  the  effectiveness  of  placing  /iPPTs  on  the  tip  mass  of  a  gravity  gradient 
boom  and  calculate  the  significance  of  their  contribution  to  meeting  system  performance 
requirements. 

The  FEM  outlined  in  Chapter  3  would  be  modified  to  incorporate  the  /iPPTs  as 
translational  inputs  located  at  the  tip  mass,  node  101.  Sensors  may  placed  at  the  tip 
mass  to  eliminate  complications  experienced  from  having  non-collocated  actuators  and 
sensors. 

An  important  consideration  in  this  research  effort  is  the  amount  of  torque  pro¬ 
duced  by  /iPPTs.  Engineering  model  testing  of  the  thrusters  measured  average  output 
values  of  25  x  1CP6  Nm.  Recall  the  limit  of  the  reaction  wheels  is  set  at  3  Nm,  a  value 
which  can  also  be  varied  by  altering  the  q/r  ratio  in  the  LQG/LTR  process.  The  design 
parameter,  R,  would  take  the  form 


which  can  be  thought  of  as 


RW 


\ 


y  mPPT 


6x3 


(5.2) 


R  = 


( 


rw  *  /3x3 


\ 


(5.3) 


y  th  *  /3X3  J 

where  rw  and  th  scaling  factors  which  can  be  modified  until  the  dynamic  response  of 
the  closed  loop  system  indicates  reaction  wheel  output  of  3  Nm  and  /rPPT  outputs  of 
25  x  10“6  Nm. 


This  analysis  could  determine  if  the  /.iPPTs  improve  or  degrade  system  perfor¬ 
mance  of  the  baseline  model.  Is  anything  gained  by  placing  sensors  on  the  tip  mass? 
Would  this  allow  improved  control  of  the  tip  mass  or  orientation  of  the  flexible  ap¬ 
pendage  relative  to  the  controlling  body?  Thrusters  used  for  station  keeping,  orbital 
maneuvers,  and  plane  changes  can  also  be  evaluated  in  a  similar  fashion  and  would 
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allow  a  combining  of  attitude  dynamics  presented  in  this  study  with  efforts  in  orbital 
dynamics. 

Appendage  dynamics  is  another  area  where  further  research  efforts  may  be  fo¬ 
cused.  The  Simulink  model  used  in  this  research  currently  plots  measured  outputs.  It 
is  possible  to  analyze  tip  mass  deflection  and  appendage  deformation  shapes  without 
measuring  the  outputs  in  Simulink.  Numerical  integrations  may  be  performed  within 
a  Matlab  program  to  propagate  the  system  forward  in  time.  The  time  history  of  the 
modal  coordinates,  77,  could  then  be  transformed  back  to  state  variables.  Plots  of  the 
physical  states  would  provide  insight  on  how  tip  mass  deflection  changes  as  the  control¬ 
ling  body  moves.  Optimal  configurations  of  the  controlling  body,  flexible  appendage, 
and  tip  mass  could  be  evaluated  against  design  and  mission  requirements  needed  for 
such  missions  where  pointing  of  the  tip  of  the  appendage  is  a  primary  concern  while  the 
controlling  body  is  only  utilized  to  reach  predetermined  pointing  requirements. 

Studying  the  motion  of  the  node  points  along  the  appendage  will  open  up  areas 
for  future  study  in  micro-meteoroid  impacts  at  various  locations  along  the  appendage. 
Also,  the  FEM  is  an  effective  testbed  in  furthering  efforts  in  technological  advances  in 
sensors  and  actuators  imbedded  within  the  appendage  material  by  providing  a  detailed 
system  dynamic  model  which  is  easily  modified  to  accommodate  placement  of  these 
sensors  and  actuators  through  the  state  space  matrices,  B  and  C. 

This  research  effort  described  in  this  document,  and  future  research  efforts  uti¬ 
lizing  the  findings  from  this  analysis,  will  strengthen  the  field  of  study  in  attitude 
control  of  small  flexible  structures.  The  near  term  utilization  of  this  study  will  support 
both  commercial  and  governmental  efforts  in  small  flexible  spacecraft  and  appendages 
constructed  from  shape  memory  composites.  Future  impacts  of  these  findings  are  as 
limitless  as  space  itself. 
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Appendix  A 


Assumed  Modes  Method 


A.l  Assumed  Modes  Model 


The  following  mathematical  development  is  based  on  the  work  of  Canavin[37] 
and  Mackison[138]  with  modifications  added  to  include  torsional  strain  storage  in  the 
flexible  appendage.  To  begin,  consider  the  satellite  system  shown  in  Figure  A.l. 


with  the  following  nomenclature: 

•  (hi,  62 ,  63):  Unit  vectors  in  the  body  frame. 

•  (01,02,03):  Unit  vectors  in  the  appendage  frame. 

•  (hj  *2, 23):  Unit  vectors  in  the  inertial  frame. 


163 


•  O:  Origin  of  the  body  frame.  Center  of  mass  of  the  undeformed  system. 

•  O':  Position  of  O  at  rest. 

•  B:  Center  of  mass  of  rigid  body. 

•  A:  Center  of  mass  of  undeformed  appendage. 

•  M *:  Total  system  mass. 

•  M:  Appendage  mass. 

•  Q:  Connection  point  of  appendage. 

•  z(t):  Motion  of  the  system  center  of  mass. 

•  R:  Vector  from  system  center  of  mass  to  beam  attachment  point. 

•  Rb ■  Vector  of  the  center  of  mass  of  the  rigid  body  with  respect  to  a  coordinate 
system  origin  in  the  rigid  body. 

•  L:  Coordinate  of  the  center  of  mass  in  the  undeformed  body. 

•  Rb  =  z  —  L 

•  Rb  =  z  —  lob  x  L 

•  ra :  Vector  from  beam  attachment  point  to  center  of  mass  of  beam. 

•  z:  Displacement  of  system  center  of  mass  from  rest. 

•  ojb'-  Inertial  angular  velocity  of  rigid  body. 

The  basis  for  the  development  includes: 

•  The  model  is  a  rigid  body  with  an  attached  flexible  cantilevered  beam. 

•  The  beam  rest  position  is  constant  relative  to  its  base. 
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•  The  motion  of  the  body  and  the  beam  consists  of  small  translations  and  rota¬ 
tions. 


No  orthogonality  requirements  have  been  placed  on  the  assumed  mode  shapes.  The 
vibration  equations  are  therefore  coupled. 

The  linearized  attitude  matrix  is 


[d]  =  [E-0] 


(A.l) 


Therefore,  the  relation  between  the  inertial  coordinates,  {i},  and  the  body  fixed  coor¬ 
dinates,  { b },  is  given  by 

b=[6]{i}  (A. 2) 


where  [0] 


( 


\ 


0  -d3  02 

[9\  =  |  #3  0  —Q\  |  (A-3) 

-e2  01  0 

The  transformation  from  the  rigid  body  frame  to  the  flexible  appendage  frame  is 


V 


J 


a  =  [ c]{b } 


(A.4) 


where  [c]  is  constant  for  an  undriven  appendage. 

For  the  undeformed  system,  the  location  of  the  center  of  mass  is  defined  by 


p  dm  =  0 


t  sys 


(A.5) 


where  p  is  the  generic  position  vector  from  the  center  of  mass  to  the  differential  mass 
element. 

Next,  the  integral  is  evaluated, 


/  p  dm  +  p  dm 

J  rigidbody  J  flexbody 


flexbody 


(A.6) 
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The  first  integral  on  the  right  hand  side  of  the  equation  is  —(M*  —  M)L  while 
the  second  of  the  integrals  is  M ( R  +  ra).  Thus, 

—  (M*  -  M)L  +  M(R  +  ra)  =  0  (A.7) 

The  inertia  dyadic  of  the  undeformed  system  is 

II*  =  IIrb  + II Ap-U  (A.8) 

where  IIrB  is  the  inertia  dyad  of  the  rigid  body  and  II%p^u  is  the  inertia  dyad  of  the 
undeformed  appendage. 

IIrb  =  Urb  +  (M*  —  M)(LLU  —  LL)  (A.9) 

H°Ap_u  =  Ilip_u  +  M[(R  +  ra)(R  +  ra)U  ~(R  +  ra)(R  +  ra)]  (A.10) 

The  kinetic  energy  for  the  system,  Tsys,  consists  of  terms  for  the  rigid  body  and 
for  the  flexible  appendage. 

Tsys  =  ^  f  v  -v  dm  +  ^  f  v  -v  dm  (A. 11) 

^  J rb  &  Japp 

where  v  is  the  inertial  velocity  of  a  generic  mass  element.  For  the  rigid  body 

TrB  =  -  f  v  -v  dm 
J  J  rb 

Trb  =  ~{M*  —  M)Rb  •  Rb  + -rob  •  Urb  '  ub  (A. 12) 

Recall 

Rb  =  z  —  L 
Rb  =  z  —  u>b  x  L 

z  =  O-O'  (A. 13) 

Inserting  Eqs.  A.  13  into  Eq.  A.  12,  the  resulting  kinetic  energy  of  the  rigid  body  becomes 

Trb  =  -(M*  -  M){z  —  ujb  x  L}  ■  {z  —  u^B  x  L}  +  -i ob  •  Urb  '  ojb 

=  ^(M*  -  M)z  ■  z  -  (M*  -  M)[z  ■  (cub  x  L)] 

“I"  7}{TI*  —  M)(ub  x  L)  ■  (u>b  x  L)  +  —lu b  ■  Urb  '  ^ b 


(A.14) 
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Using  the  Parallel  Axis  Theorem  relationships [224] 

1 1  cm  =  IIB  —  (M*  —  M)(L2I  —  L  ■  L) 

Hb  =  Hcm  +  (M*  -M)(L2I-L-L) 

IIrb  =  Urb  +  (M*  —  M)(L2I  —  L  ■  L)  (A.15) 


the  kinetic  energy  of  the  rigid  body,  with  the  inertial  dyadic  referred  to  the  center  of 
mass,  is 

Trb  =  \{M*  -  M)z  -z  +  Kjb-  II°rb  ■  coB  ~  (M*  -  M)[z  •  (lob  x  L)\  (A.  16) 


The  kinetic  energy  of  the  appendage  is 

1 


1 


Tapp  =  —  I  v  ■  v  dm  =  -  /  Tim  •  Rrn  dm 


J  app 


Japp 


where 


and 


Rm  —  z-\-R-\-ra-\-u 


Rm  —  Z  +  U  +  U>B  X  (R  +  T a) 


(A.17) 


(A.18) 


(A.19) 


Now,  expand  Eq.  A.17  and  use  77°  ,  the  dyadic  of  the  undeformed  appendage  about 
the  center  of  mass. 


T, 


1 


app 


Rm  '  Rm  dm 


Japp 


=  —  f  [z  +  u  +  U)B  X  (7?  +  ra)]  ■  [z  +  u  +  UJB  x  (R  +  ra)]  dm 

Japp 

=  -Mi  •  z  +  z  •  f  u  dm  +  i  •  lob  x  ( M(R  +  ra))  +  77/  • '«  dm 

J  Japp  J  Japp 

+  X  [  [uB  X  (7?  +  ra)]  •  [wB  X  (7?  +  ra)]  dm 

J  Japp 

+  f  u-  [lob  X  (7?  +  ra)]  din  (A. 20) 

J  app 

Consider  the  following  portion  of  Eq.  A. 20. 

X  f  [ub  X  (7?  +  ra)]  •  [wB  x  (7?  +  ra)]  dm 
J  J  aw 


(A.21) 
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Knowing  the  vector  property  a  x  b  =  —b  x  a,  and  the  following  vector  analysis [222]: 

F\  ■  (F2  x  F:i)  =  (F\  x  F2)  •  F3 
Fi  x  (F-2  x  F3)  =  (Fi  •  F3)F2  -  (F\  ■  F2)F3 

Then, 

(a  x  6)  •  (a  x  6)  =  (a  x  6  x  a)  •  b 

=  [(a  •  a)b  —  (a  •  b)a\  ■  b 

=  (a  •  a)(b  ■  b)  —  (a  •  b)(a  ■  b ) 

=  a2b 2  —  (a  •  b)2 


and  Eq.  A. 21  becomes 


\  [  uj2b(R  +  ra)2  -  {ujB  ■  ( R  +  ra))2  dm 

^  J  app 

]-u)B  (  (R  +  raj2I  ~  (R  +  ra)(R  +  r a)  dm  u>B 

^  J  aw 


(A. 22) 


(A. 23) 


The  term  in  the  integral  is  the  inertial  dyadic  of  the  undeformed  appendage  about 
the  system  mass  center,  O  (77^  u).  The  kinetic  energy  of  the  undeformed  appendage 


is 


Tav- 


p—u 


=  -Mz  ■  z  +  z  ■  f  u  dm  +  z  ■  uB  x  (A7 (7?  +  ra ))  H —  /  u  ■  ii  dm 
2  J  app  2  J  app 

+  f  u  ■  [ujb  x  (R  +  r a)]  dm  +  ]-uB  ■  II°Ap-u  '  w B  (A.24) 

The  total  system  kinetic  energy,  including  the  rigid  body  and  the  flexible  ap¬ 
pendage,  is 


'-sys 


\{M*  -  M)z  -z  +  \uB-  II°RB  •  uB  -  (M*  -  M)[z  ■  (ujb  x  L)\  +  \mz  ■  z 


+  z  ■  f  u  dm  +  z  ■  lub  x  ( M(R  +  ra))  +  -  f  u-  u  dm 

Japp  J  Japp 

+  f  u  ■  [ujb  x  (R  +  ra)\  dm  +  -uB  ■  IIAp-u  ■  uB 
J  aw  J 


(A. 25) 
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Yaw  motion 


£ 


y(r,t) 


Figure  A. 2:  Illustration  of  Spacecraft  Yaw  Producing  Torsional  Torque 


With  Eq.  A. 25,  combine  II^b  and  H%p_u  into  II* ,  combine  the  i  •  i  terms,  and 
the  moments  about  the  center  of  mass  relationship  from  Eq.  A. 7. 


'-sys 


Mz  ■  z  +  -cob  •  II*  ■  wb  +  -  f  u  ■  u  dm  +  z  ■  f  u  dm 

^  ^  J  app  J  app 


+  l ob  ■  R  x  /  u  dm  +  wg  •  /  ( ra  x  u)  dm 

J app  J app 


(A. 26) 


Neglecting  external  conservative  torques  such  as  gravity  gradient  affects,  the  po¬ 


tential  energy  of  the  system  is  due  to  the  energy  stored  in  the  deformation  of  the 
appendage.  Any  elementary  text  on  mechanics  of  materials  calculates  the  strain  energy 
of  a  deflected  beam  as[180] 


Vb’nd=\L„EI{^dr  <a-27> 

The  potential  energy  is  for  the  case  of  an  Euler-Bernoulli  beam  which  is  in  deformation 
without  any  torsional  concerns.  The  flexible  appendage  undergoing  torsion  needs  to  be 
included  in  the  equations  of  motion  of  the  system.  This  is  a  concern  since  yaw  control 
torques  will  produce  torsional  moments  along  the  longitudinal  axis  of  the  appendage 
(see  Figure  A. 2).  Roll/pitch  attitude  maneuvers  will  generate  bending  torques  and  are 
included  in  Eq.  A. 26  and  Eq.  A. 27. 

Let  7 (r,  t)  denote  the  angular  displacement  of  the  appendage.  The  angle  of  twist 
corresponding  to  a  differential  element  of  appendage  of  length  dr  is  ]dr .  Assuming 
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the  material  properties  are  uniform  in  the  appendage  and  the  angle  of  twist  is  propor¬ 
tional  to  the  torque,  the  potential  energy  of  the  appendage  in  torsion  (whose  ends  are 
not  supported  by  torsional  springs  capable  of  storing  potential  energy)  is  written  as 

Vtorsion  =  \  [  GJ(r)[^p^]2  dr  (A. 28) 

A  J  app  C/T 

where  GJ(r)  is  the  torsional  rigidity  and  Eq.  A. 28  has  the  same  structure  as  the  potential 
energy  of  a  rod  in  longitudinal  vibration.  In  addition,  if  Ipmi(r )  is  the  mass  polar  moment 
of  inertia  per  unit  length,  then  the  kinetic  energy  is  simply 

Ttorsron  =  \  f  W0[ dr  (A. 29) 

A  J  app  O* 

A. 2  Lagrangian  Equations  of  Motion 


To  determine  the  system’s  equations  of  motion,  the  Lagrangian  is  found  by  sub¬ 
tracting  the  potential  energy  components  from  the  kinetic  energy  terms  as  follows 


L  =  T  —  V 


Inserting  Eqs.  A.26-A.29  into  Eq.  A. 30,  the  system  Lagrange  becomes 


(A. 30) 


L 


u  dm 


di  <’••*)«* 


11  If  f 

=  -M  z  ■  z  + -lob  ■  II  ■  cob  +  -  /  ii-udm  +  z-  /  u 

A  A  A  Japp  Japp 

+  u>B  •  R  x  f  iidm  +  UB-  f  r  x  u  dm  +  -  f  Ipmi(r)[ 

Japp  Japp  A  Japp  Ot 

~  \S  EI02dr-\I  (A.3D 

a  Japp  Gfi  z  Japp  is i 

where,  referencing  Figure  A. 3  ,  u  is  the  flexural  displacement  and  r  is  the  coordinate 
along  the  appendage  length  to  replace  ra.  Also,  assume  the  undeformed  appendage  is 
fixed  relative  to  the  base,  ojg  =  and  small  angle  rotations,  wg  ~  {9},  II*  «  I*. 

The  Lagrangian  can  also  be  written  in  matrix  form  by  using  the  following  repre¬ 
sentations: 


2  =  {i}T{z} 


Figure  A. 3:  Displacement  of  the  Flexible  Appendage 


lob 

u 

II* 

R 

r 

R 


r 


{b}T{uB} 

wtm 

{b}TI*{b} 

{b}T{R} 

WTM 


0 

-Rs 

r2 

r3 

0 

-Ri 

-R-2 

R\ 

0 

0 

\ 

~r3 

r2 

r3 

0 

-n 

- r2 

n 

0  ) 

0 

-03 

e2  N 

03 

0 

-Oi 

-e2 

Oi 

0  ) 

and  the  basis  transformation  relationships 


m  =  mi} 
w  =  im 


which  yields 


L 


\m*{z}T{z}  +  \{6}T I* {6}  +  1  [  {«}T{u}  dm  +  {z}T[c\  I  {ft}  dm 

^  4  ^  J  app  J  app 


1 


+  {6}  [c\R  /  {«}  dm  +  {9}  [c]  /  f{u]  dm  +  -  /  Ipma  dr 

J  app  J  app  ^  J  app 


-  -  f  EI{un}T  {uu}  dr--  [  GJ{yf 

£  Japp  &  Japp 

Now,  introduce  the  distributed  coordinates 


dr 


(A. 32) 


u{r,  t)  =  J2  4>l{r)rf(t)  (A.33) 

2—1 

where  n  is  the  number  of  modes  used  to  represent  the  displacement  and 


M  =  [4>]{n} 


(A. 34) 


where  [(f)]  is  a  3xn  matrix  with  each  column  corresponding  to  a  mode  shape,  the  assumed 
mode  shapes  are  the  spatial  solutions  of  the  Euler-Bernoulli  partial  differential  equation, 
and  {77}  contains  n  modal  coordinates. 

A  similar  representation  is  done  for  the  torsional  component  where 


7  =  MM 


(A. 35) 


and  ['(/’]  is  a  lxn  torsional  mode  shape  matrix  using  the  same  modal  coordinates,  {77} 
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Substitute  Eq.  A. 34  and  Eq.  A. 35  into  Eq.  A. 32 


l  =  W{z}t{z)  +  \{<>}Tr{o} + 1  [ 

^  £  £  Japp 

+  {z} T[c]  /  {[^{v}}  dm  +  {6}T[c]R  [  {[4>\{i)}}  dm 

J app  J app 

+  {0}T[c\  f  r{[(t)\{ri}}  dm  +  ]-  [  Ipmi{[^]{v}}T {[^]{v}}  dr 

J  app  zZ  J  app 

~  \  [  EI{[<i)"]{v}}T {[(t>"]{ri}}  dr 

£  Japp 

~  If  GJ{[i’f}{ij}}T{[i’f]{r]}}  dr  (A.36) 

J*  Japp 


L  =  ^M*{z}T{z}  +  Ue}TI*{9}  +  ^{ri}T  f  [4>]T[4>\  dm{rj} 

+  {z}T[c\  [  [<j>\  dm{fj}  +  {9}t [c\R  f  [0]  dm{q} 

J  app  J  app 

+  mT[c]  f  r[</>\  dm{f]}  +  \{^}T  f  IpmA^]TW\dr{n} 

J  app  zZ  J  app 

~  \{v}T  [  EI\(j>//]T\(j)//\  dr{??} 

J*  J  app 

~  \{n}T  [  GJ[ipr]T[iJ>r]  dr{r)}  (A.37) 

J*  Japp 


x'  =  J 

'  [4>]  dm  (3xn) 

app 

(A. 38) 

X2  -  J 

f  [0]T[0]  dm  (nxn) 

app 

(A. 39) 

x3  =  J 

f  r[4>]  dm  (3xn) 

app 

(A. 40) 

Xt  -  J 

f  EI^ii]1  [0//]  dr  (nxn) 

app 

(A.41) 

xs  =  J 

f  dr  (scalar) 

app 

(A. 42) 

xt  =  J 

f  GJ[0]t[^/]  dr  (scalar) 

app 

(A. 43) 

The  Lagrangian 

can  now  be  expressed  as 

L 

\_M’{z}T{i]  +  \{6)tI'{0}  +  1{>}}TX2{>)}  +  {i}T 

[c\xi{v} 

+  {9}t[c](RX1  +  X3){r,}  +  ±{r,}TX 5{r,} 
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-  \{V}TX,{V}  ~  ^{v}TX6{r]} 


(A. 44) 


The  Lagrangian  contains  six  coordinates  which  describe  the  translation  and  ro¬ 
tation  of  the  undeformed  system  (rigid  body  motion)  and  n  modal  coordinates  which 
describe  the  twist  and  displacement  from  rest  of  the  flexible  appendage  relative  to  the 
rigid  base.  The  equations  of  motion  take  the  form[79] 


ddL  dL 
dt  dqi  dqi 


Qi 


(A. 45) 


where  Qi  are  the  generalized  forces  and  the  generalized  variables  are: 
translation:  z  (3  position  coordinates) 
rotation:  0  (3  rotational  coordinates) 
vibration:  i)  (n  modal  coordinates) 


The  partial  derivatives  of  the  Lagrangian  are: 


dL 

dz 

dL 

dz 

dL 

~dQ 

dL 

~d6 

dL 

dr) 

dL 

dr) 


0 

M*{z}  + 

0 

rw  +  HCRAi+AaH?)} 

-A4M-A6{r?} 

X2{r)}  +  X\ T[c]T{z}  +  [LiX\  +  X3)t[c}t{9}  +  X5{r,} 


The  time  derivatives  of  the  generalized  momenta  are: 


(A. 46) 
(A. 47) 
(A. 48) 
(A. 49) 
(A. 50) 
(A.51) 


ddL 

dt  dz 

=  M*{z}  +  [c\  Xx{i)} 

(A. 52) 

ddL 

~dVW 

=  I*{d}  +  [c](RXi  +  X3){rj} 

(A. 53) 

ddL 
dt  dr) 

=  X  2{f)}  +  A  f[c]T{-z}  +  (RX\  +  X3)T[c]T  {9}  +  X  5{i)} 

(A. 54) 

The  system  equations  of  motion  are  written  as 


Translation:  M*{z}  +  [c]X1{ij}  =  {Qtrans} 


(A. 55) 
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Rotation:  I*  {9}  +  [c}(RX  1  +  X3){ij}  =  {Qrot}  (A.56) 

Vibration:  V2{r/}+V5{r/}+Vf  [c]T{z}+(RX1+X3)T[c]T  {0}+X4{r]}+XQ{r]}  =  { Qdamp } 

(A.57) 

where  {Qtrans}  are  the  generalized  translational  forces  acting  on  the  system,  {Qrot} 
are  the  generalized  rotational  forces,  and  {Qdamp}  are  the  generalized  damping  forces 
within  the  flexible  appendage.  An  example  of  a  translational  system  force  would  be 
thrusters  providing  station  keeping  forces  while  attitude  control  torques,  such  as  those 
provided  by  off-axial  thrusters  or  reaction  wheels,  are  contained  within  the  rotational 
forces  category. 

The  generalized  damping  forces  within  the  flexible  appendage  will  take  the  form 
of  {Qdamp}  =  — c{?)} ,  where  the  damping  coefficient,  c,  is  not  to  be  confused  with  the 
transformation  matrix,  [c],  going  from  the  body  fixed  frame  to  the  flexible  appendage 
frame.  In  fact,  if  the  flexible  appendage  is  not  off-nominal,  meaning  the  principal  axes  of 
the  appendage  line  up  with  those  of  the  rigid  body,  then  the  transformation  matrix,  [c], 
can  be  assumed  to  equal  the  identity  matrix  and  may  be  eliminated  from  the  equations 
of  motion. 

Another  point  to  consider  is  the  dimensions  of  X2,  X4,  X5,  and  Xq.  Looking  at 
Eq.  A.57,  it  would  seem  intuitive  to  combine  these  values  together.  However,  A2  and 
X4  are  (nxn)  matrices  while  X5  and  Xq  are  scalars.  The  scalar  values  can  be  multiplied 
by  the  identity  matrix  to  put  them  in  a  form  which  will  allow  the  collection  of  similar 
terms  in  Eq.  A.57. 

Assuming  station  keeping  is  not  a  key  factor  in  the  attitude  control  conceptual 
operations  of  the  spacecraft,  and  applying  the  above  relationships,  the  equations  of 
motion  are  rewritten  as  follows: 


M*{z}  +  X1  {77}  =  0 

(A.58) 

I*  {9}  +  (RX\  +  X3){rj}  =  {Qrot} 

(A. 59) 

175 


(X2  +  Xb){f,}  +  c{i)}  +  (X4  +  X6){V}  =  -Xj {z}  -  (RX1  +  X3)T{9]  (A. 60) 

Now,  Eq.  A. 58  is  rewritten  as 

W  =  ~x: ,{,')}  (A. 61) 

and  is  substituted  into  Eq.  A. 60  to  eliminate  translational  effects  from  the  vibrational 
equation.  This  results  in  a  coupled  set  of  rotation-vibration  equations. 

p{9}  +  (Rx1  +  x3){ii}  =  {Qrot} 

(X2  +  X5  -±-XfX1){ij}  +  c{rl}  +  (XA  +  X6){ri}  =  -{RXx  +  X3)T {6\ A.62) 

The  theoretical  derivation  can  be  furthered  by  noticing  the  rotation/vibration 
coupling  terms  in  each  of  the  two  equations  are  the  transpose  of  each  other.  Let 

P  =  (RX\  +  X3)  (A. 63) 

where  P  has  dimensions  (3xn). 

Inserting  Eq.  A. 63  into  Eq.  A. 59  and  solving  for  9  yields 

9  =  (I*)-1  {Qrot}  -  (I*)~1P {fj}  (AM) 

Now,  insert  Eq.  A. 64  into  the  right  hand  side  of  Eq.  A.62. 

=  -pTm 

=  -pT((n~1{Qrot}  -  (n-'pm 
=  -PT(I*)~1{Qrot}  +  PT(P)~lP{v}  (A. 65) 

to  finally  arrive  at  the  equation  of  motion  describing  how  the  flexible  appendage  responds 
to  attitude  control  torques. 

(X‘2  +  x5-  -^Xfx,  -  PT(P)-lP){f!}  +  C{fj}  +  (X4  +  A6){??}  =  -PT(I*y1{Qrot} 

(AM) 
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where  I*  is  the  nonsingular  system  inertia  matrix  and  C  is  the  (nxn)  damping  matrix. 
Note,  the  equation  of  motion  can  be  represented  as  the  classical  dynamic  model 


M{f,}  +  C{i 7}  +  K{V}  =  F 


(A. 67) 


M  =  (X2  +  X5-j^X?X1-Pt(I*)-1P) 

K  =  (X4  +  X6) 

F  =  -PT(I*r1{Qrot}  (A. 68) 


Appendix  B 


Gravity  Gradient  Stabilization 


The  center  of  gravity  is  often  times  not  in  the  same  location  as  the  center  of  mass 
for  a  large  object  orbiting  about  a  planet.  The  differing  locations  will  result  in  a  torque 
applied  to  the  object  as  a  result  of  differing  gravitational  strength.  A  simple  solution  to 
overcome  this  applied  torque  in  a  gravity  gradient  field  is  to  attach  a  cable  to  the  satel¬ 
lite.  A  mass  is  attached  to  the  end  of  the  cable  pointing  either  toward  or  away  from  the 
earth.  The  cable  is  slowly  let  out  until  a  stable  configuration  is  achieved.  The  concept 
of  gravity  gradient  control  is  discussed  in  several  sources[80],  [90],  [159],  [176],  [238], 
and  [240]. 

B.l  Gravitational  Field 

Newton’s  Law  of  Universal  Gravitation  states: 


The  force  of  gravity  between  two  bodies  is  directly  propor¬ 
tional  to  the  product  of  their  two  masses  and  inversely  pro¬ 
portional  to  the  square  of  the  distance  between  them[186]. 


Thus,  the  gravitational  potential  energy  (V)  is  found  with 

_  GMm 


(B.l) 


where  G  =  6.6726x10  11  m3 /kgsec2 ,  M  is  the  mass  of  the  primary  body,  m  is  the  mass 
of  the  orbiting  body,  and  R  is  the  distance  separating  the  two  masses.  If  the  satellite  is 
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in  an  Earth  orbit,  then  M  =  5.9737x10 Mkg  and  the  gravitational  parameter  for  Earth 
becomes  /j  =  GM  =  3.9860xl014m3/sec2.  Substituting  the  gravitational  parameter  of 
Earth  into  Eq.  B.l  yields 


=  fxm 

R 


(B.2) 


The  gravitational  force  field  is  the  gradient  of  the  gravitational  potential  and  the 
gravitational  acceleration  acting  on  mass  m  is  defined  as: 


CLn  - 


m 


(B.3) 


Figure  B.l:  Gravity  Gradient  Torques  on  a  Near-Earth  Satellite[240] 


In  a  cartesian  coordinate  frame,  the  components  of  the  gravitational  field  are 

r«/3  =  Vq.v/3 — 
m 
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r  a/3  = 


1 

m 


d2V 

d2V 

d2V 

de'f 

de\de-2 

de\des 

d2V 

d2V 

d2V 

de\de2 

de2dez 

d2V 

d2V 

d2V 

deides 

de2de3 

Ttef" 

(B.5) 


J 


Tan  is  the  gravitational  acceleration  gradient  causing  an  acceleration  in  the  a 
direction  on  an  object  displaced  in  the  /3  direction  (a,/3  =  6i,62,6 3).  The  acceleration 
in  the  61  direction  is 

ai  =  rnei  +  ri2e2  +  Tises  (B.6) 

If  the  body  fixed  reference  frame  is  located  at  the  center  of  mass  of  the  satellite 
while  a  vector,  Ra,  points  from  the  origin  of  the  primary  body  to  the  center  of  gravity 
of  the  satellite,  the  gravitational  acceleration  in  an  inertial  frame  can  be  expressed  as 
follows: 


CLn  - 


GM 5} 

Rl 

{ R  +  r) 


R$ 


{R  +  r) 


[{R  +  r)-  {R  +  r)]  2 

-) 

R‘ 


,  ri  ,  .  ...  2 R  ■  r,  3 

os  ^t1  +  (  p)  + - — 1  2 


R 3 


R2 


Therefore,  if 


and 


then 


fJj  —*  R  •  T  — * 

-r^[-R  +  r  —  3  — j  77]  +  higher  or  derter  ms 
R6  lr 


7?  —  7?&3 


r  =  ri6i  +  r262  +  r3&3 


4  ~  +  r2^2  +  (77  -  2r3)63] 


(B.7) 


(B.8) 


(B.9) 


(B.10) 


Substituting  Eq.  B.10  into  Eq.  B.3  results  in  the  following: 

h 


rn  =  - 


R3 


(B.ll) 


180 


T22 

r33 


r 3 

2/i 

& 


and  putting  the  gradient  of  the  gravitational  field  into  matrix  form  yields 


B.2  Inertial  Gradient  Field 


(B.12) 

(B.13) 

(B.14) 


If  the  reference  frame  is  rotating  with  an  angular  velocity  cn,  then  an  inertial 
acceleration  field  is  developed.  The  gradients  from  the  inertial  field  are  added  to  the 
gravitational  gradients  to  obtain  the  total  gradient. 

The  general  expression  for  the  apparent  acceleration  of  a  point  in  a  rotating  frame 
is 

r  =  a  —  [ao  +  Cux(Cuxr)  +  Cuxr  +  2Cuxr]  (B.15) 

Only  terms  involving  r  have  gradients  which  can  be  added  to  the  gravitational 
gradients  in  Eq.  B.14.  Let 

A  =  —  [LJx(Loxr)  +  Cuxr] 

=  u)2r  —  uj(Co  ■  r)  +  rxCu 

=  A\bi  +  A2b2  +  A3b3  (B.16) 

where 

Ai  =  r2io3  -  r3io 2  +  ( ul  +  cuf  )ri  -  uiuj2r2  -  uq u3r3  (B.17) 

A2  =  r3io i  -  rio;3  +  (wf  +  u\)r2  -  w2w3r3  -  u2uiri  (B.18) 

^43  =  ri<h2  -  r2ioi  +  (ixf  +  co2)r3  -  u)3u)\ri  -  oj3ui2r2  (B.19) 


and 


to  —  ivibi  +  io2b2  +  to3b3 


(B.20) 
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After  taking  the  partial  derivatives  of  A  and  substituting  them  into  Eq.  B.14,  the 
gradient  matrix  becomes 


/ 


Ga/3  — 


\ 


(B.21) 


(rn  +  io\  +  lUj>)  ((J3  —  UJ1UJ2)  —  (cJ2  d*  uqiU3) 

—  (0)3  +  u’iu’2)  (r22  +  Lof  +  W3)  (uii  —  LO2  (-O3) 

^  (cJ2-wiu;3)  —(uii  +  w2a;3)  (T33  +  d- J 

If  a  satellite  is  orbiting  the  Earth  on  a  circular  path,  then  its  angular  velocity  is 


..2  =  JL 
0  i?3 


(B.22) 


If  the  body  reference  frame  is  orientated  such  that  b\  is  aligned  in  the  orbit  normal 
direction,  then  uj\  =  u>0  and  0J2  =  uj 3  =  0.  Inserting  these  angular  velocity  components 
into  Eq.  B.21,  and  recalling  Eq.  B.14,  the  gradient  field  for  a  circular  orbit  becomes 


^  -1  0  0  ^ 


Goi/3  ~  U0 


\ 


(B.23) 


J 


0  0  0 

0  0  3 

The  above  equation  shows  that  for  a  satellite  with  a  flexible  boom  deployed,  it 
will  experience  a  compressive  load  in  the  orbit  normal  direction  while  in  tension  along 
the  local  vertical  direction.  Thus,  a  deployed  boom  will  produce  a  stable  configuration 
either  when  the  boom  is  zenith  or  nadir  pointing. 


B.3  Gravity  Gradient  Torque 


A  gravitational  gradient  torque  acting  on  an  orbiting  body  is  expressed  as 


T ^  =  J  ( rxdg)dm 


(B.24) 


Inserting  Eq.  B.7  into  Eq.  B.24  yields 


fig)  ps  J  (rx^f-R)drn 


(B.25) 
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Letting  R  =  RE\ ,  where  E\  is  the  unit  vector  pointing  in  the  radial  direction 
from  the  primary  body,  Eq.  B.25  is  integrated  as  follows: 


7%) 


3^i  f  (R  •  r)fxR 

R 3  J  R 2 


rrdmxEi 


-j^E\x  J ( Er 2  —  rr)dm  ■  E\ 

^EhxT-Eh 


(B.26) 


where  E  =  E\E\  +  E2E2  +  E3E3  is  a  unit  dyadic  and  the  inertia  dyadic  about  the 
body’s  center  of  mass  is 


I  = 


—  rr)dm 


(B.27) 


The  gravity  gradient  with  respect  to  the  body  reference  frame  ba(a  =  1,2,3) 
becomes 


where 


f (9)  =  KE{xI  ■  Ei 


—  R  aa\ baxlapbabp  •  (iry  1  bry 


(B.28) 


Ei 


a-alba  —  an&l  +  021&2  +  (131&3 


I 

K 


lafibctbpi.Q'i  ft  1  ■  2.  3) 

3/x 

R? 


and  aa\  are  the  direction  cosines  between  the  E\  and  ba  unit  vectors. 

The  body  components  of  torque  are  found  when  the  scalar  form  of  Eq.  B.28  is 
written  as 


T(x9)  =  f(9) -6a(A  =  1,2,3) 
KaaKifnlxpe^x 


(B.29) 
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where  the  three  dimensional  permutation  tensor  is 

£a7A  =  ( baxb, 7)  •  b\(a,  P,  7,  A  =  1,  2, 3)  (B.30) 

For  the  principal  body  axes  I7/g  =  0  for  7  7^  /?.  The  torque  components  become 

t}9)  =  K(I33  -  I22)a21a31 
=  K(In  —  I33)a\ia3i 

T3(s)  =  K(I22-Iu)ana2i  (B.31) 

where  the  aa/g  terms  are  still  the  direction  cosines  (i.e.  ba  =  aapEp). 


B.4  Equations  of  Motion 


Euler’s  equation  of  motion  for  a  satellite  in  a  circular  orbit  with  orbital  angular 
velocity,  lus,  is 

h  +  CS8xh  =  T{ [g)  (B.32) 

where 

h  =  u>il\bi  +  u>2I2b2  +  ui3I3b3 
f(g)  =  + 


Three  sequential  rotations  0\  about  the  b\  axis,  02  about  the  b-2  axis,  and  03 
about  the  b3  axis  are  used  to  describe  the  orientation  of  the  satellite’s  principal  axes 
with  respect  to  the  orbiting  reference  frame.  The  body  reference  frame  is  expressed  in 
terms  of  the  fixed  reference  frame  as  follows: 


\ 


b2 


=  R(63)R(62)R(e  1) 


E2 


IV 


(B.33) 
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where 


R(Oi)  = 


1  0  0 
0  c6\  s9\ 

0  —  s6\  c6\ 


(B.34) 


R(e2)  = 


c6 2  0  -S0-2 


0  1 


s9 2  0  C0-2 


(B.35) 


and 


R(03)  = 


R123  =  R(93)R(02)R(9i) 

( 


\ 


C0 3  S0 3  0 

-S0  3  C0  3  0 


0 


0  1 


(B.36) 


C02C03 

C03S01S02  +  C0\S03 

—  c9\C03S02  +  S0\S03 

-C02S03 

—  S01S02S03  +  C01C03 

C0\S02S03  +  S0\C03 

S0  2 

—  C02S01 

C0\C02 

(B.37) 


) 


The  expanded  form  of  the  body  components  of  the  satellite  angular  velocity,  us, 


is 


coi  =  9i+  ojo{c93s9is02  +  c9\s93) 

U2  =  02+  UO(c9iC03  -  S01S02S03 ) 

U>3  =  03  +  LUO(—S0iC02)  (B.38) 

with  loq  =  -p-.  Linearizing  yields 


—  9\+  ^003 

L02 

=  02  +  Wo 

U>3 

=  03  —  WQ01 

(B.39) 
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The  component  equations  in  scalar  form  are 

=  I\LOi  +  —  h) 

T. 2^  =  I 2^2  +  W\Uz(I\  —  I3) 

T3(9)  =  I^  +  u^h-h)  (B.40) 


Taking  the  derivative  of  Eq.  B.39  and  substituting  into  Eq.  B.40  produces 

T A  =  h0i  +  U0O3)  +  02  +  ^o)03  ~  woOi)(h  ~  h) 

=  I2O2  +  01  +  ^0^3)03  —  W0&l)(Il  —  I3 ) 
r3(9)  =  Iz(ez-uQOl)  +  0l+Uoe3)02  +  uo){l2-h)  (B.41) 

Recall  the  body  components  of  the  gravity  gradient  torque,  t[9\  T^9\  and  T39\ 
shown  in  Eq.  B.31.  Using  small  angle  approximations,  the  torque  components  are  now 
expressed  as 


-t(s) 


-1  (s) 


vu 


(9) 


3(Un 


0 


(/1  -  h)6 2 


(B.42) 


)  \  (A  ~  h)63  j 

Inserting  these  torques  into  Eqs.  B.41,  the  linearized  equations  of  motion  for  small 


angular  deviations  become 


0  —  h0i  +  W(A)  +  (I2  —  ~  ^(A) 

(B.43) 

0  =  I282  +  3ca02(/3  —  /i)#2 

(B.44) 

0  =  h03  ~  cucA)  +  {h  —  U)(4:W q03  +  ^cA) 

(B.45) 

In  addition,  the  gravity  gradient  restoring  torques  resulting  from  small  angular 
deviations  are  obtained  by  neglecting  the  small  8  coupling  terms. 


T\  =  -u%(I2-I3)d  1 
T2  =  —  3Wq(/3  —  Ii)02 
T3  =  -4wl{h-h)8z 


(B.46) 
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note,  the  restoring  torques  vanish  for  a  symmetrical  satellite. 

B.5  Stability  Considerations 

The  derivations  provided  in  the  previous  section  were  for  a  yaw  (9 1),  pitch  (#2), 
roll  (63)  configuration.  However,  FalconSAT-3  uses  a  roll  (#1),  pitch  (62),  yaw  (63) 
convention.  To  avoid  confusion  between  research  results  presented  in  this  document  and 
those  found  by  the  FalconSAT-3  design  teams,  Eqs.  B.43-B.45,  the  linearized  equations 
of  motion  of  a  rigid  body  in  a  circular  orbit  (roll,  pitch,  and  yaw,  respectively)  are  as 
follows: 


h(9 1  —  <^0^3)  +  (h  —  ^i)(4wq^i  +  <^0^3)  =  0 

(B.47) 

1 20  2  +  3cJo(/l  —  13)02  =  0 

(B.48) 

h(@3  +  wo^i)  +  (h  —  4i)(wq^3  —  wo^i)  =  0 

(B.49) 

Because  the  pitch-axis  equation,  Eq.  B.48,  is  decoupled  from  the  roll/yaw  equa¬ 
tions,  Eq.  B.47  and  Eq.  B.49,  (i.e.  there  are  no  9\  or  63  terms)  consider  the  characteristic 
equation  of  the  pitch  axis  given  by 

s2  +  3co%^-— —  =  0  (B.50) 

h 

If  1 1  >  I3,  then  the  characteristic  roots  are  pure  imaginary  numbers  and  the  pitch 
equation  is  a  simple  harmonic  oscillator.  If  I\  <  I3,  then  one  of  the  characteristic 
roots  is  a  positive  real  number  and  the  pitch  equation  is  unstable.  For  an  unstable 
configuration,  the  satellite  will  swing  away  from  the  equilibrium  configuration  when 
disturbed.  Therefore,  it  is  required  that  I\  >  I3  for  pitch  stability. 

For  roll/yaw  stability  analysis,  Eq.  B.47  and  Eq.  B.49  are  rewritten  as 

9\  +  (k\  —  \)ujq03  +  4u)qA;i$i  =  0 

9'i  +  (1  —  &3)wo$i  +  ^0^363  =  0 


(B.51) 
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where 


h 

k2 


h~h 

h 

h~h 

h 


Eqs.  B.51  can  be  combined  into  one  coupled  roll/yaw  system  equation  as 


f1  °) 

(V 

0  uj0(k  1  -  1) 

U] 

4wq/ci 

0  \ 

+ 

+ 

v°  V 

V  h  ) 

^  w0(l  -  k3)  0 

V  63  ) 

V  0 

j 

v  / 

(B.52) 


The  characteristic  equation  of  Eq.  B.52  becomes 


s2  +  SLOo(ki  —  1) 

s<^o(l  —  k3)  s2  +  u%k3 


(B.53) 


or 


s4  +  (1  +  3Aq  +  kik3)uQS2  +  4kik3ujQ  =  0 


(B.54) 


The  roots  of  the  coupled  roll/yaw  system  characteristic  equation  become  pure 
imaginary  numbers  if  and  only  if 


k\  k3  >  0 
1  +  3 k\  +  k\  k3  >  0 
(1  +  3/ci  +  k±k3)2  —  \d>k\k3  >  0 


(B.55) 


The  requirements  in  Eqs.  B.55,  combined  with  the  pitch  stability  requirement  of  I\  >  I3, 
govern  the  complete  stability  of  the  gravity  gradient  satellite  equilibria.  This  stability 
result  is  illustrated  using  a  stability  diagram  in  the  {k\ ,  k3)  plane  as  shown  in  Figure  B.2. 
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Figure  B.2:  Gravity  Gradient  Spacecraft  Stability  Regions 


Appendix  C 


Rotational  Dynamics 


C.l  Rotational  Kinematics 

When  studying  the  equations  of  motion  of  an  object  one  must  understand  both 
kinetics  and  kinematics.  Kinetics  is  the  study  of  the  forces  acting  on  a  body  while 
kinematics  describes  the  motion  of  that  body.  This  section  will  concern  itself  with 
the  theory  of  rotational  kinematics.  Kinetics  will  be  discussed  in  Section  C.2.  ’’The 
subject  of  rotational  kinematics  is  somewhat  mathematical  in  nature  because  it  does  not 
involve  any  forces  associated  with  motion”  [238].  The  three  mathematical  approaches 
to  identifying  rotational  motion  are  direction  cosines,  Euler  angles,  and  quaternions. 

C.1.1  Direction  Cosines 

The  first  mathematical  method  of  describing  rotational  kinematics  is  direction 
cosines.  In  a  rigid  body  there  exists  a  body  reference  frame  A.  This  frame  consists  of 
a  right-hand  set  of  three  orthogonal  unit  vectors  {01,02,03}.  Another  reference  frame 
B  will  contain  a  different  set  of  three  right-handed  orthogonal  unit  vectors  {61 , 62 ,  ^3}  - 
The  dot  product  law  is  used  to  describe  one  reference  frame  in  terms  of  the  other. 

a  ■  b  =  \a\\b\  cos  lab  (C.l) 

Using  the  dot  product  B  can  now  be  expressed  in  terms  of  A  as  follows: 

b\  =  Clio}  +  6*1202  +  6*1303 
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t>2  ~  C2 1  U 1  +  C*22®2  +  C23O3 

^3  =  ^3101  +  C 32^2  +  <^3303 


(C.2) 


where  C\j  =  bi  ■  dj  is  the  cosine  of  the  angle  between  bi  and  dj  and  Cij  is  called  the 
direction  cosine. 

An  equivalent  way  of  expressing  the  directional  cosines  is  in  matrix  form. 


i  r  \ 


h 

b2 

\h  j 


( 


\  l  .  \ 

a  1 
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) 


a2 

ya3  J 


(C.3) 


Cn  C 12  C13 
C21  C22  C23 
C31  C32  C33 

The  direction  cosine  matrix  is  also  called  the  rotation  matrix  or  the  coordinate  trans¬ 
formation  matrix  [2]. 

The  kinematical  equations  of  Poisson  describe  the  functional  relationships  of  the 
direction  cosines  and  their  rates  [44], 


C\  1  =  C12W3  —  C 13W2 
Cl  2  =  C13UJ1  —  CuUJ3 
C\3  =  C11U2  ~  Ci2iO\ 

C‘2\  =  C22W3  —  C23CU2 
C22  =  C23u)i  —  C21CU3 

C*23  =  C21OJ2  —  C*22^1 

C'i\  =  C32W3  —  C33U2 
C32  =  C33U11  —  C31U3 

C33  =  C31LO2  —  C32UJ1  (C.4) 


The  direction  cosines  are  easy  to  calculate  but  do  require  integrating  9  equations  to 
solve  rotational  kinematics.  Also,  the  direction  cosines  are  not  intuitive  since  the  values 
are  not  expressed  as  angles. 
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C.1.2  Euler  Angles 

A  more  intuitive  way  of  looking  at  rotational  kinematics  is  with  Euler  angles. 
These  angles  are  easier  to  visualize  than  direction  cosines.  A  perfect  example  of  Euler 
angles  is  the  yaw,  pitch,  and  roll  of  an  aircraft.  Euler  angles  involve  rotating  about  the 
three  axes  of  the  body  as  shown  in  Figure  C.l. 


e3 


For  each  rotation  about  an  axis,  a  rotational  matrix  is  calculated.  Consider  a 
rotation  of  a  body  axis  e  with  respect  to  a  reference  frame  E  as  shown  in  Figure  C.2. 
The  components  of  E  along  the  e  directions  are  calculated  as 


&i  =  Ei 

e-2  =  E2  cos  0\  +  E3  sin  6\ 

=  —E2  sin  9\  +  E3  cos  6\ 


63 


(C.5) 
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or  in  matrix  form  as 


where  c6  and  s9  represent  cos#  and  sin0. 


(C.6) 


Figure  C.2:  Sequential  orthogonal  rotations  of  the  e  reference  frame  about  the  E  refer¬ 
ence  frame  [44] 


A  rotation  about  each  axis  can  be  represented  as 


(C.T) 


(C.8) 
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=  R(93) 
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V  "3  ) 

\  ^  J 

where 


(C.9) 
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0  c6\  s6\ 
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-s03  c#3  0 


(C.10) 


(C.ll) 


(C.12) 


V  0  0  V 

Referring  again  to  Figure  C.l,  assume  the  body  is  undergoing  a  rotation  first 
about  the  e\  axis  followed  by  a  rotation  about  e2  and  then  e3.  The  body  reference 
frame  is  expressed  in  terms  of  the  fixed  reference  frame  as  follows: 


or 


(  -  ^ 

(  *  \ 

ei 

62 

=  R(03)R(02)R(9  i) 

E2 

V  "3  ) 

V  E:i  / 

(C.13) 


#i23  =  R(e3)R{d2)R(d  i) 
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c92c9s 

c93s9is9-2  +  c9\s93 

—c9ic93s92  +  s9\s93 

-c9-2s93 

— s9is62s93  +  c9\c93 

c9\s92s93  +  s9\c93 

s92 

—c92s9i 

C#i C$2 

The  angular  velocities  for  this  rotation  sequence  are 

i0\  =  0  cos  (p  +  ip  sin  9  sin  (p 

L02  =  ip  sin  9  cos  cp  —  9  sin  (p 
L03  =  (j)  +  ip  cos  9 

Solving  for  <p,  9,  ip  yields 


(C.14) 


(C.15) 


d>  ~  —  7i  —  sin<bcos#  —  cos  0  cos  0  sin  9  t oo  (C.16) 

^  9  y  y  cos  (p  sin  9  —  sin  cp  sin  9  0  J  y  ui3  ^ 

These  three  equations  appear  easy  to  integrate  but  trigonometric  functions  require 
greater  computational  time  than  addition  and  multiplication.  Although  Euler  angles 
have  fewer  equations  to  integrate  than  direction  cosines,  a  larger  amount  of  processing 
time  is  used  to  calculate  the  angular  velocities.  The  only  concern  is  the  singularity  in 
Eq.  C.16  when  9  equals  0°.  For  a  gravity  gradient  stabilized  satellite,  this  singularity 
becomes  a  concern.  To  avoid  the  singularity  in  pitch,  a  roll-pitch-yaw  (1-3-2)  rotation 
sequence  is  used.  This  rotation  sequence  yields 


ip  cos  9  —  cos  ip  sin  9  sin  ip  sin  9  u\ 

cp  =  -7,  0  cos  ip  —sin  ip  ijjo  (C-17) 

^  cos  9  ' 

^  9  y  y  0  sin  ip  cos  9  cos  ip  cos  9  y  y  013 

Now  the  singularity  occurs  when  the  pitch  angle  equals  90°  and  gravity  gradient  dy¬ 
namics  of  the  satellite  will  prevent  this  singularity  from  occuring. 
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C.1.3  Quaternions 

The  third  mathematical  approach  to  identifying  rotational  motion  is  using  quater¬ 
nions.  Instead  of  the  trigonometric  functions  of  Euler  angles,  quaternions  use  algebraic 
relations.  This  provides  quicker  computations  and  also  eliminates  the  possibility  of  a 
singularity  appearing.  An  Euler  axis,  also  known  as  a  principal  axis,  is  an  axis  that 
is  fixed  in  the  body  frame  and  is  stationary  in  the  inertial  frame.  This  axis  is  special 
because  any  combination  of  rigid  body  rotation  is  described  as  a  single  rotation  about 
the  Euler  axis.  The  eigenaxis  vector  e  =  (ei,e2,e3)  is  simply  the  direction  cosines  of 
the  Euler  axis  relative  to  both  the  body  and  inertial  reference  frames. 

The  quaternions,  Euler  parameters,  are  then  defined  as 

91  =  ei  sin(<9/2) 

92  =  e2sin(0/2) 

93  =  e3  sin(<9/2) 

94  =  cos(0/2)  (C.18) 


where  0  is  the  rotation  angle  about  the  Euler  axis.  After  applying  trigonometric  iden¬ 
tities  the  rotation  matrix  becomes 


/ 


\ 


1  -  2(<?2  +  <?i)  %1<?2  +  9394)  2(9193  -  q2q4) 

R  =  I  2(^1^2-9394)  1- 2^  +  93)  2(^293  +  9194)  (C.19) 

^  2(^3  +  q2q4)  2(q2q3  -  9^4)  1  -  2 (qf  +  qf)  ) 

Substituting  Eq.  C.19  into  Eqs.  C.4  and  solving  for  angular  velocity  we  obtain 


ui  =  2(9194  +  9293  -  9392  -  949l) 

u>2  =  2(9294  +  9391  -  9193  -  9492) 

W3  =  2(9394  +  9192  -  9291  -  9493) 


(C.20) 
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Applying  a  constraint  equation  such  that 

1i  +  I2  +  li  +  94  =  1  (C.21) 

and  differentiating  it 

0  =  2(qiqi  +  q2q-2  +  1.3 1.3  +  1414)  (C.22) 

We  can  now  combine  Eqs.  C.20  and  Eq.  C.22  into  matrix  form  as  follows: 
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Rearranging  Eq.  C.23  results  in  the  kinematic  differential  equations  for  quater¬ 


nions  as 
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12 
13 

V  I4  / 


(C.25) 


— ^3  0  u)  1  lo2 

lo2  — tcq  0  W3 

— ui\  2  —  C03  0 

As  a  result  of  the  fewer  number  of  equations  to  integrate  compared  to  direction 

cosines,  and  the  absence  of  singularities  and  trigonometric  functions,  modern  spacecraft 
orientation  is  now  commonly  described  in  terms  of  quaternions  [238]. 


C.2  Rigid  Body  Dynamics 

Kinematic  differential  equations  are  useful  for  integrating  the  motion  of  a  satellite. 
An  understanding  of  kinetics  is  required  to  study  how  this  satellite  reacts  to  forces  acting 
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upon  it.  This  section  will  first  discuss  the  properties  of  angular  momentum,  moments 
of  inertia,  and  how  moments  will  affect  the  equations  of  motion.  The  section  concludes 
with  a  discussion  on  stability  requirements  for  both  linear  and  nonlinear  systems. 

C.2.1  Kinetics 

Just  as  a  moving  body  has  translational  momentum  comprised  of  its  velocity  and 
mass,  a  rotating  body  will  also  have  rotational  momentum.  This  rotational  momentum 
is  referred  to  as  angular  momentum.  Similar  to  its  counterpart,  angular  momentum  is 
comprised  of  angular  velocity,  oj,  and  inertia,  I.  Thus  the  total  angular  momentum  is 
expressed  as 

H  =  10  (C.26) 

Both  H  and  Co  are  3x1  vectors  so  I  is  a  3  x  3  matrix  known  as  the  inertia  matrix. 

A  rigid  body  does  not  necessarily  have  to  possess  uniform  density  distribution. 
It’s  density  can  vary  with  position  relative  to  the  system’s  center  of  mass.  The  moment 
of  inertia  matrix  contains  all  of  the  information  of  the  mass  distribution  within  a  rigid 
body  [240]. 

/ ( y 2  +  z2)dm  —  f  xydm  —  f  xzdm 

—  f  xydm  f(x2  +  z2)dm  —  J  yzdm  (C.27) 

—  /  xzdm  —  J  yzdm  J (x2  +  y2)dm  j 

This  inertia  matrix  only  needs  to  be  recalculated  if  the  mass  distribution  alters  from  its 
original  configuration. 

The  kinetic  relation  between  moments  and  angular  momentum  is 

M  =  H  (C.28) 

Using  a  body  fixed  reference  frame  Eq.  C.28  becomes 

M  =  —H  +  ujx  H 
dt 


(C.29) 
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Applying  the  matrix  form  of  the  cross  product  yields 


M  =  Ilo  +  loxIlo 

where  w1  is  a  scew-symmetric  matrix  of  the  form 
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(C.30) 
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To  further  simplify  the  equations  of  motion,  assume  a  principal  axis  frame.  The 
inertia  matrix  becomes  a  diagonal  matrix 
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and  Eq.  C.30  becomes 
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which  reduces  to  Euler’s  rotational  equations  of  motion  for  a  rigid  body. 

Mi  =  JiT’i  —  ( J2  —  J3)uj  2^3 
M2  =  J2^>2  —  {J3  —  Jl)^iU>3 
M3  =  J3lo3  —  (Ji  —  J2)u\u2 


(C.34) 


These  three  equations  are  coupled,  nonlinear  ordinary  differential  equations.  The 
rotational  motion  of  a  rigid  body  can  be  completely  described  when  Euler’s  rotational 
equations  of  motion  are  combined  with  the  kinematic  differential  equations  from  Sec¬ 


tion  C.l. 
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C.2.2  Stability 

Stability  is  a  major  concern  with  satellite  control  theory.  An  unstable  satellite 
may  drift  away  from  its  operational  orbital  path,  disrupt  communications,  or  waste 
excessive  fuel  and  reduce  its  mission  effectiveness.  An  object  near  an  equilibrium  point 
is  considered  stable  when  small  disturbances  result  in  small  changes.  To  understand 
stability  about  an  equilibrium  point  consider  a  pendulum.  Two  equilibrium  points 
exist.  One  point  occurs  when  the  pendulum  is  pointing  straight  down  and  the  other 
point  when  the  pendulum  is  pointing  straight  up.  At  both  equilibrium  points  the 
pendulum  will  remain  motionless  if  no  external  forces  are  applied.  However,  only  one 
of  these  points  are  stable.  A  small  disturbance  applied  to  the  pendulum  while  pointing 
down  will  cause  it  to  oscillate  slightly  about  its  equilibrium  point.  This  configuration  is 
considered  stable.  A  small  disturbance  applied  to  the  pendulum  while  it  is  pointing  up 
will  cause  a  large  change  in  its  orientation.  This  point  is  unstable. 

Current  theory  supports  two  concepts  of  stability.  The  first  concept  deals  with 
Lagrange  stability.  Lagrange  stability  applies  when  small  disturbances  result  in  bounded 
changes.  Liapunov  stability  occurs  when  an  object  resting  at  an  equilibrium  point 
experiences  small  disturbances  and  only  small  changes  are  detected.  Now  consider  the 
pendulum  example  previously  mentioned.  When  the  pendulum  is  pointing  downward  it 
is  considered  to  be  both  Liapunov  and  Lagrange  stable.  When  the  pendulum  is  pointing 
straight  up,  any  small  change  will  result  in  a  drastic  change  in  location.  Thus,  this  point 
is  Liapunov  unstable.  But  the  angle  of  the  pendulum  is  bounded  (i.e.  it  may  oscillate 
slightly  about  the  equilibrium  position)  so  the  system  is  considered  Lagrange  stable. 
From  this  it  is  apparent  that  a  Liapunov  stable  system  is  also  Lagrange  stable  but  a 
Lagrange  stable  system  is  not  necessarily  Liapunov  stable. 

To  test  if  a  linear  system  is  considered  stable  a  concept  from  control  theory  is 
applied.  ”A  linear  system  is  said  to  be  stable  if  and  only  if  the  roots  of  its  characteristic 
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equation  have  negative  real  parts”  [44].  This  is  the  same  as  saying  that  the  real  part  of 
the  eigenvalues  of  the  system  are  only  negative.  A  zero  root  is  considered  asymptotically 
unstable.  A  system  is  asymptotically  stable  if  changes  from  an  equilibrium  point  result 
in  returning  to  the  same  equilibrium  point. 

Liapunov  developed  two  theorems  to  determine  if  a  nonlinear  system  is  stable. 
The  first  theorem  states  that  a  nonlinear  system  is  asymptotically  stable  in  a  region 
near  an  equilibrium  point  if  and  only  if  its  linearized  approximation  is  asymptotically 
stable.  The  second  theorem  states  if  a  positive  definite  function  can  be  found  such  that 
its  derivative  is  negative  definite  the  nonlinear  system  is  stable  [2].  This  function  is 
often  referred  to  as  the  Liapunov  function.  For  a  more  indepth  discussion  of  stability 
and  methods  to  determining  a  Liapunov  function  refer  to  [44],  pages  113-126. 


Appendix  D 


FalconSAT-3 


At  the  same  time  CTD  was  conducting  tests  on  the  CoilAble  boom,  the  United 
States  Air  Force  Academy  (USAFA)  was  beginning  the  design  phase  for  their  fourth 
satellite,  FalconSAT-3.  The  Academy’s  Space  Systems  Research  Center  (SSRC)  sup¬ 
ports  the  Department  of  Defense  (DoD)  Defense  Planning  Guide  mandate  to  ’’establish 
and  sustain  a  cadre  of  space  professionals”  by  providing  an  opportunity  for  AFA  cadets 
to  ’’learn  space  by  doing  space”  within  a  safe,  supervised  environment.  This  allows  the 
cadets  to  gain  real-world,  practical  experience  working  with  small  satellite  and  launch 
systems.  Students  in  the  fields  of  astronautics,  physics,  math,  computer  science,  and 
management  are  supervised  by  members  of  AFA  faculty  and  contractor  support  in  the 
research  areas  of  systems  engineering,  sounding  rocket  systems,  and  micro/nanosat  sys¬ 
tems. 

The  FalconSAT-3  spacecraft  is  a  0.46m  cube  weighing  approximately  50  kg.  A 
solid  model  drawing  of  the  current  spacecraft  is  shown  in  Figure  D.l.  The  spacecraft  is 
manifested  for  launch  aboard  an  Atlas  V  Evolved  Expendable  Launch  Vehicle  (EELV) 
in  October,  2006[207].  While  deployed  at  an  altitude  of  560knr  and  an  inclination  of 
35.4°,  the  spacecraft  is  designed  to  provide  an  on-orbit  platform  that  supports  three 
DoD  experiments:  the  Flat  Plasma  Spectrometer  (FLAPS),  the  Plasma  Local  Anom¬ 
alous  Noise  Environment  (PLANE),  and  the  Micro  Propulsion  Attitude  Control  System 
(MPACS). 
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Figure  D.l:  A  Solid  Model  Drawing  of  FalconSAT-3 
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The  objective  of  the  FLAPS  experiment  is  to  characterize  the  effects  of  non- 
Maxwellian  charged  particle  distributions  on  formation,  propagation,  and  decay  of 
ionospheric  plasma  bubbles[208].  Successful  completion  of  this  objective  will  support 
the  validation  of  the  plasma  bubble  and  radio  wave  scintillation  nowcasting  and  fore¬ 
casting  system  associated  with  DoD’s  Communication/Navigation  Outage  Forecasting 
System  (C/NOFS).  The  FLAPS  experiment  will  make  fine  resolution  (10%)  measure¬ 
ments  of  ionospheric  plasma  spectra  differential  in  energy  and  angle  and  could  support 
simultaneous  multi-point  in  situ  measurements  of  a  single  plasma  bubble  structure  with 
another  satellite.  FLAPS  requires  FalconSAT-3  to  provide  ±5°  attitude  control  in  the 
ram  direction  (2-axis)  with  attitude  knowledge  of  ±1°  for  post  processing  of  the  data. 

The  PLANE  experiment’s  objective  is  to  characterize  the  plasma  turbulence  in 
the  space  environment  surrounding  a  satellite [209].  PLANE  is  designed  to  distinguish 
the  turbulence  in  the  ambient  environment  (global  effects)  from  variations  in  the  plasma 
population  co-moving  with  the  satellite  (local  effects).  PLANE  is  also  capable  of  quan¬ 
tifying  plasma  perturbations  caused  by  active  systems  on  the  satellite,  such  as  a  bring 
propulsion  system.  The  experiment  does  not  require  state  vector  data  for  real  time 
operations.  However,  data  post  processing  requirements  call  for  ±20°  pointing  in  the 
ram  direction  (2-axis),  a  rate  of  change  in  the  pitch  and  yaw  axis  of  less  than  2  deg/sec, 
and  attitude  knowledge  of  ±3°. 

The  MPACS  thruster  is  based  on  the  micro  pulsed  plasma  thruster  (rnPPT) 
developed  at  Air  Force  Research  Lab’s  Propulsion  Directorate  over  the  last  several 
years[210].  A  three-electrode  version  of  the  mPPT  was  developed  by  Busek  for  advanced 
engineering  development  on  the  FalconSAT-3  flight  (see  Figure  D.2).  Flight  heritage  of 
the  Micro  Propulsion  Attitude  Control  System  (MPACS)  thruster  is  to  be  established 
on  FalconSAT-3  in  the  hopes  of  demonstrating  that  the  thruster  can  be  used  as  an 
effective  actuator  as  part  of  a  microsatellite  ADCS  [211].  The  experiment’s  objective  is 
to  operate  MPACS  as  part  of  the  FalconSAT-3  ADCS  for  100  cumulative  hours.  Data 
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Figure  D.2:  3-Axis  MPACS  Cluster 


concerning  the  beginning  and  end  of  life  performance  characteristics,  as  well  as  detection 
of  any  adverse  operational  interactions  between  MPACS  and  the  satellite,  need  to  be 
quantified.  To  demonstrate  attitude  stabilization  functions,  four  3-axis  clusters  are 
mounted  on  the  nadir  plate  of  the  satellite  (Figure  D.3)  and  one  2-axis  cluster  is  placed 
on  the  tip  mass  of  the  gravity  gradient  boom  (Figure  D.4)[205]. 

The  attitude  determination  requirements  imposed  by  the  MPACS  experiment  are 
nebulous  at  best.  The  calculated  moments  of  inertia  of  FalconSAT-3  with  an  opera¬ 
tionally  deployed  boom  are  shown  in  Table  D.l.  The  FalconSAT-3  ADCS  design  team 
are  conducting  simulation  efforts  to  determine  the  impact  the  approximately  100  fi/N 
of  thrust  produced  by  MPACS  will  have  on  effectively  controlling  the  satellite. 

FalconSAT-3  is  the  USAFA’s  first  attempt  at  designing  a  3-axis  stabilized  and 
controlled  spacecraft.  Aware  of  the  power  and  volume  limitations  inherent  with  small 
satellites,  the  ADCS  team  is  utilizing  both  active  and  passive  actuators  in  the  design. 
Attitude  sensors  on  the  spacecraft  consist  of  four  Space  Quest  model  SS-256  (2-axis) 
sun  sensors  and  one  Billingsley  TFM100G2  Fluxgate  Magnetometer.  Three  Space  Quest 
magnatorque  rods  (5.0  Am2)  provide  active  control  while  CTD  is  developing  an  EMC 
gravity  gradient  boom  for  passive  control.  The  EMC  longerons  will  be  both  the  pri¬ 
mary  deployment  mechanism  and  the  principal  structural  members  of  the  boom.  Al- 
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Figure  D.3:  Four  MPACS  Clusters  on  the  Nadir  Plate  of  FalconSAT-3 


Figure  D.4:  A  2-Axis  MPACS  Cluster  Attached  to  the  EMC  Boom 
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Table  D.l:  Moments  of  Inertia  with  Deployed  Boom  Configuration 


MOI 

kg-rn2 

lbm-in2 

lbf-ft-s2 

lbf-in-s2 

I XX 

67.40 

23.0£’4 

49.72 

596.60 

lyy 

67.45 

23.1774 

49.74 

596.98 

Izz 

1.31 

4.468773 

0.96 

11.57 

though  the  application  of  composite  materials  has  the  potential  to  reduce  weight  and 
improve  performance,  the  cost  is  experiencing  increased  deflection  and  vibration  within 
the  structure [241].  With  the  flexible  characteristics  of  the  boom  (1.5  Hz  bending  fre¬ 
quency  and  1.7  Hz  torsional  frequency)  there  is  great  concern  of  the  impact  this  flexible 
structure  will  have  on  the  control  of  a  small  satellite. 


